Large Scale Network Balance Under Stress and Real-life Affective Stress Reactivity

Introduction

The following code is specific to the task related fMRI task analysis of the Stress Resilience and the Brain in Medical Students (STRAIN-MD) study. Details of the study can be found in an earlier publication. In short, the study consisted of a within subject design examining real-life stress exposure with Ecological Momentary Assessments(EMA) measures, and an fMRI stress exposure using the Socially Evaluated Cold Pressor task(SECPT). The main aim of the study as whole is to link changes is large scale brain networks under stress to real-life stress and resilience measures. Details of the real-life stress procedure can be found in the previously mentioned publication. During the fMRI session, participants engaged in three tasks. The fMRI tasks are designed to activate specific parts of the brain. Primarily, the Executive Control (ECN), Default Mode (DMN), and Salience (SN) networks. The ECN task is a 2-Back where participants respond every time they see 2-Back number. The DMN task is an associative retrieval task. For this task, participants were shown images that moved to 1 of four corners of the screen outside the scanner. They were then asked to indicate where they saw the images using button presses in the scanner. The third and final task that targeted the SN was an emotional oddball. For this task, participants had to respond to a stream of faces whenever a novel face with an emotional expression was presented. In addition to the actual tasks and their contrasts, we also collect peripheral physiological data. This includes

• Heart Rate
• Skin Conductance
• Salivary Cortisol

These measures were used to validate our laboratory stress induction procedure. In this notebook, we validate the stress induction procedures using the above heart rate and saliva data. We additionally derive a lab-based stress reactivity measure from the salivary cortisol response to stress during the SECPT. We also derive a measure from the real-life stress data that reflects stress reactivity/resilience in daily life. Both the lab stress measure, and the real-life stress measure can then be used in subsequent analysis of fMRI related activity and overall task performance in the corresponding sections below. In order to replicate these analyses, the below libraries are use, in addition to a file with custom functions supplied in the repository.

library(knitr)
knitr::opts_chunk$set(include=TRUE, results='show', message=FALSE, warning=FALSE, comment=NA, echo=T, tidy=TRUE, concordance=TRUE, autodep=TRUE)
options(knitr.table.format = "html", booktabs=TRUE) 
#Cleaning Libs
library(readxl)       # Read excel sheets
library(tidyverse)  
library(plyr)
library(dplyr)
library(broom)         # tidy functions
library (hablar)
library(janitor)
library(data.table)
library(zoo)
# Basic Stats
library (psych)       # Nice descriptives and scaling
library (psycho)      # We use the dPrime Tool
library(Hmisc)        # Girl I dont event know
library(corrr)        # Correlation
library(ppcor)        # Partial Correlation
library(corrplot)     # Correlation matrix plotting
#Linear Modeling
library(lmerTest)     # Main Mixed Models
library(afex)         # Advanced Options Mixed Models
library(optimx)       # Optimizer for MMs
library(jtools)       # Uses summ function for readable models
library(performance)  # Model Diagnostics
library(emmeans)      # used for follow-up tests on lmers
#Plotting
library(ggplot2)      # Standard plotting
library(ggpubr)       # Arrange plots
library(sjPlot)       # Plot regressions and lmer
library(car)          # I use the qqPlot function here
# Fonts and Asthetics
library(stargazer)
library(extrafont)
library(kableExtra)
# Functions and cores
library(future)       # Parallel Processing
library(foreach)      # Parallel for each functins
library(parallel)
library(doSNOW)       # do function over cores
source("functions.R") # Custom Functions
# Make sure lmer contrasts are set to sum 0
options(contrasts=c("contr.sum", "contr.poly"))
# Assign Cores for Parallel processing
ncores <- as.numeric(availableCores())-1
# Load fonts
#font_import()
loadfonts(device="pdf")
# Themes used for ggplots
ggtheme1 <- theme (text=element_text(size=16),
                  legend.position = "none",
                  axis.text.y = element_blank(),axis.text.x = element_text(size=16),
                  axis.ticks.y = element_blank(),
                  plot.title = element_text(size=16, hjust=0.5),
                  panel.background = element_rect(fill="transparent"),
                  panel.grid.minor.y = element_line(size=3),
                  panel.grid.major = element_line(colour = "aliceblue"))
ggtheme2 <- theme(text=element_text(size=11, family="ArialMT"),
                  axis.line = element_line(size = 1, colour = "grey"),
                  panel.background = element_rect(fill="transparent"),
                  plot.background = element_rect(fill = "transparent", color = NA), # bg of the plot
                  panel.grid.minor = element_line(colour="grey95"),
                  panel.grid.major = element_line(colour = "grey95")) 
fmri_colors <- c("#00BFC4", "#7CAE00", "#F8766D", "#00B8E7", "#0CB702", "#E68613" ) # DMN, ECN, SN, Early then Late
stress_colors <-  c("#00BFC4",  "#F8766D", "#7CAE00","#00B8E7", "#E68613")               

1 Real-life Stress

We first start with the stress induction in real life, since to some extent this is a replication of previously done analyses. Here, we used an exam and control week as stressor to see its impact on mood an physiology. This was done prior to any scanning, and you can find out more details in the publication. In short, we found that our stress weeks result in increased subjective stress, and that affect and physiology both change in relation to stress exposure, and subjective stress reports. Our primary goal for the sake of the current project is to derive a measure of stress reactivity in real life that we can use as a fixed effect in our fMRI analyses. We should be able to use this measure to link effects of real life to the lab, to fMRI task performance, or to fMRI task related activity and most importantly large scale network balance.

To this end, we do an interpretation of the residual based regressor, with subjective stress as our exposure measures, and positive affect as our resilience component. This is based on a few different studies (Rafael Kalisch’s work, and some from Sanne Boij). We contemplated using something like emotional inertia and the instability, but these measures are only useful for between subject designs, rather than within.

1.1 Data

We load in the file that we derive in our paper mentioned above. This consists of processed EMA items for subject stress during a week with an exam and another week without an exam (IE. Stress and Control weeks respectively). All data cleaning has already been performed, so we need minimal processing at this stage. The full details of the processing pipeline that was used to generate this data can be found in our previous study, with the associated notebook made available on GitHub. All we need to do is load the data, and do some simple retyping.

# Import Data = Retype
df_EMA <- read.csv2("data/EMA_Clean_Anon.csv", sep = ",")
df_EMA <- df_EMA %>% retype() #%>% print()
df_EMA$week_type <- df_EMA$week_type
df_EMA$Week_Type <- factor(df_EMA$week_type, levels=c(0,1), labels=c("Control", "Stress"))

We also want to dichtomoize our subjective stress levels. We do this by assessing the overall subjective stress, and determining subject-level means. If the stress is higher that the mean, then we get a stress situation. If it is lower, then we get a no stress situation.

df_EMA <- df_EMA %>% dplyr::group_by(castor_record_id) %>% dplyr::mutate(ss_tot=(event_tot+activity_tot+social_tot), ss_tot_c=scale(ss_tot,center=TRUE), 
                                                ss_tot_m=mean(ss_tot_c, na.rm=TRUE), ss_tot_bin=factor(ifelse(ss_tot_c>ss_tot_m, "Stress", "Non-Stress")) ) %>% ungroup() 
df_EMA$mean_stress_z <- scale(df_EMA$mean_stress, center=TRUE, scale=T)

We additionally want to come up with a measure from these EMA weeks that can be related to neuroticism or depressive personality to make sure that whatever measure we derive is actually valid and related to an already established measure of stress sensitivity.

df_psych <- read.csv2("/project/3013068.02/stats/EMA/Questionnaires.csv", sep=",") %>% dplyr::rename(sub_id=castor_record_id)

1.2 Statistics

Primary Models: Daily life differences

We first need to get the main models where we look at the differences in subjective stress (SS) and positive affect (PA) between the control and stress weeks. We use mixed effects models, with random slopes and intercepts for each subject. This allows us to actually get subject level fixed effects for the week type. Its similar to doing individual regressions for each subject to get the Betas for each one. These Betas are used in the next stage. Note that the results of these model are counter intuitive because the subject level fixed effects we for week type is set as Control-Stress. We therefore need to reverse the signs of the fixed effects in order to actually get the results we want (Stress-Control).

Subjective Stress

Model

contrasts(df_EMA$Week_Type) <- contr.treatment(c("Control", "Stress"))
model.stressBIN_week <- glmer(ss_tot_bin ~ Week_Type + (1+Week_Type|castor_record_id), contrasts=list(Week_Type="contr.treatment"),
                              data=df_EMA, family=binomial(), control=glmerControl(calc.derivs=FALSE))
asis_output(tab_model(model.stressBIN_week, show.stat=TRUE, show.se=TRUE)$knitr)

	ss tot bin
Predictors	Odds Ratios	std. Error	CI	Statistic	p
(Intercept)	0.74	0.04	0.67 – 0.82	-5.64	<0.001
Week_Type [Stress]	1.67	0.14	1.42 – 1.98	6.02	<0.001
Random Effects
σ2	3.29
τ00 castor_record_id	0.11
τ11 castor_record_id.Week_TypeStress	0.38
ρ01 castor_record_id	-1.00
N castor_record_id	83
Observations	6153
Marginal R2 / Conditional R2	0.020 / NA

Plot

df_EMA.ss_sum <- df_EMA %>% dplyr::group_by(castor_record_id, week_type, ss_tot_bin) %>% filter(ss_tot_bin=="Stress") %>% dplyr::summarise(ss_tot_bin=length(ss_tot_bin))

`summarise()` has grouped output by 'castor_record_id', 'week_type'. You can override using the `.groups` argument.

plot.stress_week <- ggplot(df_EMA.ss_sum, aes(x=week_type, y=ss_tot_bin, color=castor_record_id)) + 
           geom_line(alpha=0.15, stat="summary") + 
            geom_smooth(aes(color=NA), method="lm", alpha=1, size=1, se=F)+
           scale_x_continuous(breaks=c(0,1),labels = c("Control\nWeek","Stress\nWeek"))+
           ylab("Number of beeps with stress") +
           ggtheme2 +
           theme(legend.position = "none", 
                 axis.title.x=element_blank(), axis.text.x=element_text(size=10, colour="black"),
                 axis.title.y=element_text(size=10), axis.text.y=element_text(size=8, colour="black" ))
plot.stress_week

Positive Affect

Model

model.pos_week <- lmer(mood_positive_cs ~ Week_Type + (1 + Week_Type|castor_record_id),
                       data=df_EMA, 
                       control=lmerControl(calc.derivs =F),
                       contrasts=list(Week_Type=contr.treatment(2)) )
asis_output(tab_model(model.pos_week, show.stat=TRUE, show.se=TRUE)$knitr)

	mood positive cs
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	6.45	0.10	6.25 – 6.64	63.73	<0.001
Week_Type2	-0.77	0.11	-1.00 – -0.55	-6.80	<0.001
Random Effects
σ2	3.67
τ00 castor_record_id	0.74
τ11 castor_record_id.Week_TypeStress	0.86
ρ01 castor_record_id	-0.44
ICC	0.17
N castor_record_id	83
Observations	6155
Marginal R2 / Conditional R2	0.033 / 0.196

Plot

plot.pa_week <- ggplot(df_EMA, aes(x=week_type, y=mood_positive_cs, color=castor_record_id)) + 
            geom_line(alpha=0.15, stat="summary") + 
           geom_smooth(aes(color=NA), method="lm", alpha=1, size=1, se=F)+
           scale_x_continuous(breaks=c(0,1),labels = c("Control\nWeek","Stress\nWeek"))+
           ylab("Average positive affect (a.u.)") +
           ggtheme2 +
           theme(legend.position = "none", 
                 axis.title.x=element_blank(), axis.text.x=element_text(size=10, colour="black"),
                 axis.title.y=element_text(size=10), axis.text.y=element_text(size=8, colour="black" ))
plot.pa_week

Reisdual Resilience Score

We now extract the fixed and random effects of each of our models, so that we can use them in a single linear model. We first do this for SS, and then PA. Finally we combine them into a single data frame that we can use for the analysis. We do this by extracting the fixed effects estimates from the main models of subjective stress and positive affect for the random effect of subject.

# Subjective stress
fixed_stress2 <- coef(model.stressBIN_week)$castor_record_id
fixed_stress2$id <- rownames(fixed_stress2)
random_stress2 <- ranef(model.stressBIN_week)$castor_record_id
# Positive affect
fixed_pos <- (coef(model.pos_week)$castor_record_id)
fixed_pos$id <- rownames(fixed_pos)
random_pos <- ranef(model.pos_week)$castor_record_id

We finally combine them into a single data frame that we will use in the secondary model

df_stressrx <- merge(fixed_stress2, fixed_pos, by="id", suffixes=c("_stress","_pa"))
names(df_stressrx) <- c("id","stress_int","stress_fe",  "pa_fe", "pa_int", "pa_mean")

Model

We now go to our second level models. Here, we model the effect of changes in subjective stress reports on changes in mood between the two weeks. We see a significant relationship between the two. The more stress you experience, the lower your positive affect gets relative to baseline. We extract from the model of change in affect as a function of change in subjective stress reports the fitted values, residuals, predicted value, and plot them. From the residuals, we see that those above the line are people who experience less change in affect in response to stress, and those below experience more reductions in positive affect in response to stress. So those with negative residuals would be stress sensitive, and those without would be resilient.

model.residual_regressorsBIN <- lm(pa_fe ~ stress_fe, data=df_stressrx)
asis_output(tab_model(model.residual_regressorsBIN, show.stat=TRUE, show.se=TRUE)$knitr)

	pa fe
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.43	0.12	0.20 – 0.66	3.70	<0.001
stress_fe	-0.84	0.16	-1.16 – -0.51	-5.11	<0.001
Observations	83
R2 / R2 adjusted	0.244 / 0.234

# Binarized Model
df_stressrx$fitsBIN <- model.residual_regressorsBIN$fitted.values
df_stressrx$resBIN <- model.residual_regressorsBIN$res
df_stressrx$effBIN <-model.residual_regressorsBIN$effects
df_stressrx$predBIN <- predict(model.residual_regressorsBIN)

Plot

plot.residual <- ggplot(df_stressrx, aes(x=stress_fe, y=pa_fe)) + 
    geom_smooth(method="lm", se=F, colour="grey")+                          # Plot linear relationship
    geom_point(aes(y=fitsBIN), shape=1, size=0.2)+                          # Fits dots on the line
    geom_segment(aes(xend=stress_fe, yend = fitsBIN), alpha = .2) +      # Fits the lines betweem the fitted values on the line, and residuals
    geom_point(aes(color=resBIN)) +                                         # Draw residual dots + colour them
    scale_color_gradient2(low = "red", mid = "ghostwhite", high = "blue")+  # Scale the colour
    xlab("FE Stress") + ylab("FE Positive Affect") +
    ggtheme2 +
    theme(legend.position = "none", axis.title.x=element_text(size=12), axis.title.y=element_text(size=12),axis.line.x=)
plot.residual

Full Plot

ggarrange(
          ggarrange(NULL,plot.stress_week,  NULL, plot.pa_week, ncol=2, nrow=2, widths=c(0.05,1,1), labels=c("A","", "B")), 
          NULL, 
          plot.residual, widths=c(1,0.05,1.2), 
          ncol=3, labels=c("","C", "")); ggsave("figures/Figure_2_EMA.svg", device="svg", dpi=300, bg="white")

No summary function supplied, defaulting to `mean_se()`
`geom_smooth()` using formula 'y ~ x'
Removed 468 rows containing non-finite values (stat_summary).No summary function supplied, defaulting to `mean_se()`
`geom_smooth()` using formula 'y ~ x'
Removed 468 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'

Correlation with Psych

Finally we check the correlation between our subjective stress measure, and our personality traits. This allows us to validate that our stress resilience score is related to previously linked scores from the literature. To do this, we merge the psych data frame containing the questionnaires that participants filled in, with the residual data frame, and do some magic to put all the data together. We then plot a correlation matrix. Our scores are linked to neuroticism, which is great!

df_psych.sub <- df_psych %>% dplyr::rename(id=sub_id) %>% dplyr::select(id, ("neo_o"), contains("neo_c"), contains("neo_e"), contains("neo_a"), contains("neo_n"), contains("bdi_total"), "trait_total")
df_stress.sub <- df_stressrx %>% dplyr::select(id, contains("fe"), starts_with("res") )
df_stress.psych <- merge(df_stress.sub, df_psych.sub, by="id")%>% retype()
names(df_stress.psych) <-c("id", "PA", "SS_BIN", "Res_BIN", "NEO_O", "NEO_C", "NEO_E", "NEO_A", "NEO_N", "BDI", "STAI_T")
df_stress.psych.cor <- df_stress.psych[c(2:11)]
df_stress.psych.cor <- data.frame(lapply(df_stress.psych.cor, function(x) scale(x, center = TRUE, scale=T)))

cormat.str <- rcorr(as.matrix(df_stress.psych.cor))
corrplot(cormat.str$r, method="number", p.mat=cormat.str$P, type="lower", diag=FALSE, insig="pch", cl.pos="n")

asis_output(tab_model(lm(Res_BIN ~ 1 + NEO_N, data=df_stress.psych ))$knitr)

	Res BIN
Predictors	Estimates	CI	p
(Intercept)	0.33	-0.05 – 0.71	0.085
NEO_N	-0.21	-0.43 – 0.00	0.054
Observations	81
R2 / R2 adjusted	0.046 / 0.034

ggplot(df_stress.psych, aes(y=NEO_N, x=Res_BIN)) + geom_smooth(method="lm") + ggtheme2

2 Laboratory Stress Induction

Before we do any analysis regarding the effects of stress on the brain and our tasks, we need to validate our stress induction procedure. Since this is a within subject design, all participants underwent both the stress induction procedure, and a matched control procedure. The stress procedure involved a novel experimenter placing the participants foot in ice water of 3 minutes, followed by a difficult mental arithmetic task for another 3 minutes. The control procedure was performed by a familiar experimenter, with room temperature water, and an easy arithmetic task. To validate the procedure, we first analyze our collected physiology data which includes saliva and heart rate measures. We also collect skin conductance data, but that still needs to be done, and is outside the scope of the current paper.

2.1 Data

The saliva measures are derived from samples collected during the scanning at 0 min, 6 min, 30 min, 50 min, and 120 min (approximately) from stress induction/control procedure. Additional samples are collected in the stress/control weeks at home. These samples are analysed offsite and returned in an excel sheet, with the timing logs being recorded in Castor during the scanning. The heart rate measures are also collected during this time via Brainamp with a pulse oxymeter. These measures were cleaned in matlab using HeRa, and then compiled using a custom script. We first do some simple cleaning procedures to make analysis easier for each of the measures.

Saliva

Salivary cortisol and amylase is collected several times in our study. Twice during an exam week, and a control week (non-exam), in addition to the 5 saliva samples during each of the MRI sessions (stress and control). A spreadsheet with the values of the saliva assays was sent from the offsite analysis site. We also record logs for the time at which each sample was acquired. Since this is our main stress validation measure, we need to make sure to get the files clean first. We therefore import the logs containing timing information for our samples to preprocess them and add them to the results from the saliva assays. A few things are done here, including subsetting the data frames, and calculating the relative time of each samples using the first sample as the starting time. We can then import the actual results from the salivary assays. One participant had an incorrect time stamp, so we also correct it here (it had the wrong hour value).

# Get castor df
df_castor <- read_excel("data/castor_data_export.xlsx",sheet = "Study results")

-
/
                                                                                                                                                                                

# Get demographics DF
df_demographics <- read_delim("data/demo_anon.csv", delim="\t") %>% clean_names() %>% dplyr::rename(subject_number="subject_id") %>%
  dplyr::filter(!is.na(contraceptive_use)) %>% 
  dplyr::mutate(sex=ifelse(contraceptive_use=="Male", "Male", "Female"),
                contraceptive_use=ifelse(contraceptive_use=="c", NA, contraceptive_use)) 

indexing demo_anon.csv [====================================================================================================================================] 401.92MB/s, eta:  0s
                                                                                                                                                                                  

                                                   
# Clean Subject numbers
df_castor <- df_castor %>% separate("Record Id", c("subject_number", "dbl_sess"), "_s")
df_castor$dbl_sess[is.na(df_castor$dbl_sess)]  <- 1
df_castor <- subset(df_castor, df_castor$dbl_sess != 2) 
# Subset df
df_castor <- subset(df_castor, select=c("subject_number", "mri_data_control", "mri_data_control_1",
                                                  "cort_1", "cort_2", "cort_3", "cort_4", "cort_5", 
                                                   "cort_1_1", "cort_2_1", "cort_3_1", "cort_4_1", "cort_5_1"))
df_castor[df_castor < 0] <- NA
df_castor <- df_castor[df_castor$mri_data_control!="01-01-2999",] 
df_castor <- unique(df_castor)
# Format dates
df_castor$mri_date_control <- as.Date(df_castor$mri_data_control, format="%d-%m-%Y")
df_castor$mri_date_stress <- as.Date(df_castor$mri_data_control_1, format="%d-%m-%Y")
# Calculate which is first scan
df_castor$first_scan_diff <- NULL
df_castor$first_scan_diff <- df_castor$mri_date_control - df_castor$mri_date_stress
df_castor$first_scan_diff <- as.numeric(df_castor$first_scan_diff, units="days")
# Depednding on difference betweem dates, set first scan
df_castor$first_scan <- "Control"
df_castor$first_scan[df_castor$first_scan_diff > 0] <- "Stress"
# Calcualte time from first sample
df_castor$t_1_c <- as.numeric(difftime(as.ITime(df_castor$cort_1), as.ITime(df_castor$cort_1), units='min'))
df_castor$t_2_c <- as.numeric(difftime(as.ITime(df_castor$cort_2), as.ITime(df_castor$cort_1), units='min'))
df_castor$t_3_c <- as.numeric(difftime(as.ITime(df_castor$cort_3), as.ITime(df_castor$cort_1), units='min')) 
df_castor$t_4_c <- as.numeric(difftime(as.ITime(df_castor$cort_4), as.ITime(df_castor$cort_1), units='min')) 
df_castor$t_5_c <- as.numeric(difftime(as.ITime(df_castor$cort_5), as.ITime(df_castor$cort_1), units='min')) 
# Stress 
df_castor$t_1_s <- as.numeric(difftime(as.ITime(df_castor$cort_1_1), as.ITime(df_castor$cort_1_1 ), units='min')) 
df_castor$t_2_s <- as.numeric(difftime(as.ITime(df_castor$cort_2_1), as.ITime(df_castor$cort_1_1 ), units='min')) 
df_castor$t_3_s <- as.numeric(difftime(as.ITime(df_castor$cort_3_1), as.ITime(df_castor$cort_1_1 ), units='min'))
df_castor$t_4_s <- as.numeric(difftime(as.ITime(df_castor$cort_4_1), as.ITime(df_castor$cort_1_1 ), units='min')) 
df_castor$t_5_s <- as.numeric(difftime(as.ITime(df_castor$cort_5_1), as.ITime(df_castor$cort_1_1 ), units='min')) 
# Make long
df_castor.long <- pivot_longer(df_castor, cols=t_1_c:t_5_s, names_to=c("sample_number"), values_to="time")
# Rename Sample Number
## Control Samples
df_castor.long$sample_number[df_castor.long$sample_number=="t_1_c"] <- "_CMR_1"
df_castor.long$sample_number[df_castor.long$sample_number=="t_2_c"] <- "_CMR_2"
df_castor.long$sample_number[df_castor.long$sample_number=="t_3_c"] <- "_CMR_3"
df_castor.long$sample_number[df_castor.long$sample_number=="t_4_c"] <- "_CMR_4"
df_castor.long$sample_number[df_castor.long$sample_number=="t_5_c"] <- "_CMR_5"
## Stress samples
df_castor.long$sample_number[df_castor.long$sample_number=="t_1_s"] <- "_SMR_1"
df_castor.long$sample_number[df_castor.long$sample_number=="t_2_s"] <- "_SMR_2"
df_castor.long$sample_number[df_castor.long$sample_number=="t_3_s"] <- "_SMR_3"
df_castor.long$sample_number[df_castor.long$sample_number=="t_4_s"] <- "_SMR_4"
df_castor.long$sample_number[df_castor.long$sample_number=="t_5_s"] <- "_SMR_5"
# Get saliva results + clean up
df_saliva <- read_excel ("data/BaLLS_CortisolSamples_Master.xlsx", sheet="mastersheet")

-
/
                                                                                                                                                                                

df_saliva <- df_saliva %>% retype() %>% clean_names()
df_saliva$subject_number <- tolower(df_saliva$subject_number)
df_saliva <- merge(df_saliva, df_demographics, by="subject_number")
# Merge with castor df
df_saliva <- merge(df_saliva, df_castor.long, by=c("subject_number","sample_number"), all=T)
df_saliva <- rename_(df_saliva, "first_scan"="first_scan.y")
#Variable for Scan Type
df_saliva <- df_saliva %>% dplyr::mutate(scan_type=ifelse(str_detect(sample_number, "CMR"), "Control", NA), scan_type=ifelse(str_detect(sample_number, "SMR"), "Stress", scan_type))
# Correct the wrong time
df_saliva$cort_3_1[df_saliva$cort_3_1=="12:25" & df_saliva$subject_number=="sub_023"] <- "13:25"

We next calculate the responses of interest from the cortisol data. These will be single values of responsivity that we can use in follow-up models and include the AUC and Henckens response. For the AUCg and AUCi, we follow the formula derived from Pruessner et al. (2003). We then measure the differences in both the AUCi and the AUCg in the stress and control sessions. For the HPA response from Henckens et al., we instead calculate the difference between the two peak measures in salivary cortisol. The expected peak in our study is sample 3, so we look at the difference in sample 3 in the stress vs control sessions, while correcting for baseline acquired at home. This baseline correction is superior because it is independent of order effects, and scanner apprehension.

# Make a wide df
df_saliva.wide <- df_saliva %>% dplyr::select(subject_number,  first_scan, sample_number, cortisol, alpha_amylase, time) %>% 
  pivot_wider(id_cols=c(subject_number), names_from=sample_number, values_from=c(cortisol, alpha_amylase, time), names_sep="", values_fn=list(cortisol=mean, alpha_amylase=mean, time=mean)) 
# Get the AUCg 
## AUCg Control
df_saliva.wide$AUCg_control <- (
    (((df_saliva.wide$cortisol_CMR_1+df_saliva.wide$cortisol_CMR_2)/2)*(df_saliva.wide$time_CMR_2)) +
(((df_saliva.wide$cortisol_CMR_2+df_saliva.wide$cortisol_CMR_3)/2)*(df_saliva.wide$time_CMR_3-df_saliva.wide$time_CMR_2)) +
(((df_saliva.wide$cortisol_CMR_3+df_saliva.wide$cortisol_CMR_4)/2)*(df_saliva.wide$time_CMR_4-df_saliva.wide$time_CMR_3))
)
## AUCg Stress
df_saliva.wide$AUCg_stress <- (
    (((df_saliva.wide$cortisol_SMR_1+df_saliva.wide$cortisol_SMR_2)/2)*(df_saliva.wide$time_SMR_2)) +
(((df_saliva.wide$cortisol_SMR_2+df_saliva.wide$cortisol_SMR_3)/2)*(df_saliva.wide$time_SMR_3-df_saliva.wide$time_SMR_2)) +
(((df_saliva.wide$cortisol_SMR_3+df_saliva.wide$cortisol_SMR_4)/2)*(df_saliva.wide$time_SMR_4-df_saliva.wide$time_SMR_3)))
# Now get AUCi
## AUCi Control
df_saliva.wide$AUCi_control <- df_saliva.wide$AUCg_control - (df_saliva.wide$cortisol_CMR_1* ((df_saliva.wide$time_CMR_2)+(df_saliva.wide$time_CMR_3-df_saliva.wide$time_CMR_2)+(df_saliva.wide$time_CMR_4-df_saliva.wide$time_CMR_3))) 
## AUCi Stress
df_saliva.wide$AUCi_stress <- df_saliva.wide$AUCg_stress - (df_saliva.wide$cortisol_SMR_1* ((df_saliva.wide$time_SMR_2)+(df_saliva.wide$time_SMR_3-df_saliva.wide$time_SMR_2)+(df_saliva.wide$time_SMR_4-df_saliva.wide$time_SMR_3))) 
## AUC Differences
df_saliva.wide$AUCg_difference <- df_saliva.wide$AUCg_stress - df_saliva.wide$AUCg_control
df_saliva.wide$AUCi_difference <- df_saliva.wide$AUCi_stress - df_saliva.wide$AUCi_control
# Henckens
df_saliva.wide$HPA_act <- (df_saliva.wide$cortisol_SMR_3 - df_saliva.wide$cortisol_CMR_1)/(rowMeans(cbind(df_saliva.wide$cortisol_CW_1, df_saliva.wide$cortisol_CW_2, df_saliva.wide$cortisol_SW_1, df_saliva.wide$cortisol_SW_2), na.rm=TRUE))

In order to plot the results, we also derive the mean times at which the samples we acquired. Finally we can also clean up some factors and other related items in the data frame to make analysis easier. Also filter out the non-scan related saliva samples because these are not really of interest here.

#Time Points for Scan Samples
df_saliva$time_point[df_saliva$sample_number=="_CMR_1"|df_saliva$sample_number=="_SMR_1"]<- mean((df_saliva.wide$time_CMR_1+df_saliva.wide$time_SMR_1)/2, na.rm=TRUE)
df_saliva$time_point[df_saliva$sample_number=="_CMR_2"|df_saliva$sample_number=="_SMR_2"]<-mean((df_saliva.wide$time_CMR_2+df_saliva.wide$time_SMR_2)/2, na.rm=TRUE)
df_saliva$time_point[df_saliva$sample_number=="_CMR_3"|df_saliva$sample_number=="_SMR_3"]<-mean((df_saliva.wide$time_CMR_3+df_saliva.wide$time_SMR_3)/2, na.rm=TRUE)
df_saliva$time_point[df_saliva$sample_number=="_CMR_4"|df_saliva$sample_number=="_SMR_4"]<-mean((df_saliva.wide$time_CMR_4+df_saliva.wide$time_SMR_4)/2, na.rm=TRUE)
df_saliva$time_point[df_saliva$sample_number=="_CMR_5"|df_saliva$sample_number=="_SMR_5"]<-mean((df_saliva.wide$time_CMR_5+df_saliva.wide$time_SMR_5)/2, na.rm=TRUE)
# Make into wide
df_saliva.wide_sub <- df_saliva.wide %>% dplyr::select(subject_number, AUCg_control:HPA_act) 
# Add to the main df
df_saliva <- merge(df_saliva, df_saliva.wide_sub, by="subject_number") 
# Factor variables and filter non needed data
df_saliva$scan_type <- factor(df_saliva$scan_type,levels=c("Control","Stress"), labels=c("Control", "Stress"))
df_saliva.scan <- df_saliva %>% filter(!is.na(first_scan))

Heart Rate

Heart rate measures were collected using a Brainamp based system. We load in the heart rate data, and start by cleaning it up. The actual denoising and artifact removal is done in Matlab using HeRa (mentioned above). For this analysis, we simply import the collected results from this program, and then proceed with a few simple cleaning steps. First we import the data. Then we clean. Certain subjects have issues with the naming of their fMRI EEG recordings of HR data (for example, due to stopping the scan and restarting, or a program crash). So we need to do some cleaning for the names and sessions for each of these subjects individually. This is done below. Its a little rudimentary but this was coded in the first days of using R and its too cumbersome to edit now. Basically, if it ain’t broke, don’t fix it. Finally, we re-order our session indicators, and double check to make sure they all have factors that let us plot and analyse the data effectively.

# Load data
df_HR <- read.csv("data/fmri_HR_clean.csv")
# Drop duplicates
df_HR<-df_HR[!duplicated(df_HR$File_Name),]
#Sub 10:
df_HR$Run[df_HR$SUB_NR == 10 & df_HR$Run == 5] <- 4
df_HR$Run[df_HR$SUB_NR == 10 & df_HR$Run == 0] <- 5
df_HR$Run_time[df_HR$SUB_NR == 10 & df_HR$Run == 4] <- 29
df_HR$Run_time[df_HR$SUB_NR == 10 & df_HR$Run == 5] <- 45
# Sub 52 s1_r1 errors:
df_HR<-df_HR[!(df_HR$SUB_NR==52 & df_HR$Run==1 & df_HR$Session==2),]
df_HR$Run[df_HR$SUB_NR == 52 & df_HR$Run == 2 & df_HR$Session==2] <- 1
df_HR$Run[df_HR$SUB_NR == 52 & df_HR$Run == 3 & df_HR$Session==2] <- 2
df_HR$Run[df_HR$SUB_NR == 52 & df_HR$Run == 4 & df_HR$Session==2] <- 3
df_HR$Run[df_HR$SUB_NR == 52 & df_HR$Run == 5 & df_HR$Session==2] <- 4
df_HR$Run[df_HR$SUB_NR == 52 & df_HR$Run == 0 & df_HR$Session==2] <- 5
# SUB67  Session 2
df_HR$Run[df_HR$File_Name=="sub67_s2_r1-2_fMRI_run_1_TR1500ms_189vols_-10-294s_hera.mat"] <- 2
df_HR$Run[df_HR$File_Name=="sub67_s2_r1-2_fMRI_run_2_TR1500ms_475vols_629-1361s_hera.mat"] <- 3
df_HR$Run[df_HR$File_Name=="sub67_s2_r1-2_fMRI_run_3_TR1500ms_200vols_1380-1700s_hera.mat"] <- 4
df_HR$Run[df_HR$File_Name=="sub67_s2_r1-2_fMRI_run_4_TR1500ms_475vols_1730-2462s_hera.mat"] <- 5
# SUB 71 Session 1
df_HR$Run[df_HR$File_Name=="sub_71_s1_r2_fMRI_run_1_TR1500ms_475vols_161-893s_hera.mat"] <- 6
df_HR$Run[df_HR$File_Name=="sub_71_s1_r2_fMRI_run_2_TR1500ms_200vols_915-1235s_hera.mat"] <- 7
df_HR$Run[df_HR$File_Name=="sub_71_s1_r2_fMRI_run_3_TR1500ms_475vols_1267-1999s_hera.mat"] <- 8
df_HR$Run[df_HR$File_Name=="sub_71_s1_r2_fMRI_run_4_TR1500ms_200vols_2019-2339s_hera.mat"] <- 9
df_HR$Run[df_HR$File_Name=="sub_71_s1_r2_fMRI_run_5_TR1500ms_475vols_2371-3103s_hera.mat"] <- 10
df_HR$Session[df_HR$Session==0 & df_HR$SUB_NR==71] <- 2
# SUB 72 Session 1
df_HR$Run[df_HR$File_Name=="sub_72_s1_r1_fMRI_run_2_TR1500ms_200vols_2900-3220s_hera.mat"]=2
df_HR$Run[df_HR$File_Name=="sub_72_s1_r1_fMRI_run_3_TR1500ms_475vols_3615-4348s_hera.mat"]=3
df_HR$Run[df_HR$File_Name=="sub_72_s1_r1_fMRI_run_4_TR1500ms_200vols_4396-4716s_hera.mat"]=4
df_HR$Run[df_HR$File_Name=="sub_72_s1_r1_fMRI_run_5_TR1500ms_475vols_4764-5497s_hera.mat"]=5
df_HR$Run[df_HR$File_Name=="sub_72_s1_r2_fMRI_run_1_TR1500ms_475vols_176-908s_hera.mat"]=6
df_HR$Run[df_HR$File_Name=="sub_72_s1_r2_fMRI_run_2_TR1500ms_200vols_944-1264s_hera.mat"] =7
df_HR$Run[df_HR$File_Name=="sub_72_s1_r2_fMRI_run_3_TR1500ms_475vols_1289-2021s_hera.mat"]=8
df_HR$Run[df_HR$File_Name=="sub_72_s1_r2_fMRI_run_4_TR1500ms_200vols_2049-2369s_hera.mat"]=9
df_HR$Run[df_HR$File_Name=="sub_72_s1_r2_fMRI_run_5_TR1500ms_475vols_2409-3142s_hera.mat"]=10
df_HR$Session[df_HR$Run>0 & df_HR$SUB_NR==72] <- 1
# SUB 72 Session 2
df_HR$Run[df_HR$File_Name=="sub_72_s2_r1_fMRI_run_1_TR1500ms_475vols_1005-1737s_hera.mat"]=1
df_HR$Run[df_HR$File_Name=="sub_72_s2_r1_fMRI_run_2_TR1500ms_200vols_1768-2088s_hera.mat"]=2
df_HR$Run[df_HR$File_Name=="sub_72_s2_r1_fMRI_run_3_TR1500ms_475vols_2421-3154s_hera.mat"]=3
df_HR$Run[df_HR$File_Name=="sub_72_s2_r1_fMRI_run_4_TR1500ms_200vols_3171-3491s_hera.mat"]=4
df_HR$Run[df_HR$File_Name=="sub_72_s2_r1_fMRI_run_5_TR1500ms_475vols_3505-4238s_hera.mat"]=5
df_HR$Run[df_HR$File_Name=="sub_72_s2_r2_fMRI_run_1_TR1500ms_475vols_374-1106s_hera.mat"]=6
df_HR$Run[df_HR$File_Name=="sub_72_s2_r2_fMRI_run_2_TR1500ms_200vols_1133-1453s_hera.mat"]=7
df_HR$Run[df_HR$File_Name=="sub_72_s2_r2_fMRI_run_3_TR1500ms_475vols_1526-2258s_hera.mat"]=8
df_HR$Run[df_HR$File_Name=="sub_72_s2_r2_fMRI_run_4_TR1500ms_200vols_2337-2657s_hera.mat"]=9
df_HR$Run[df_HR$File_Name=="sub_72_s2_r2_fMRI_run_5_TR1500ms_475vols_2686-3418s_hera.mat"]=10
df_HR$Session[df_HR$Session==0 & df_HR$SUB_NR==72] <- 2
# Now remove anything extra that wasnt cleaned up:
df_HR$Run[df_HR$Run==0] <- NA
# Factor to indicate late and early
df_HR$time_point <- 1 
df_HR$time_point[df_HR$Run > 5] <- 2
df_HR$time_point <- factor(df_HR$time_point, levels=c('1','2'),   labels=c('Early','Late'))
# Factir fir variable for stress and control
df_HR$Session <- factor(df_HR$Session, levels=c('1','2'), labels=c('Control','Stress'))
# Lets also add the first scan coloumn from castor:
df_HR_First <- df_castor
df_HR_First <- separate(df_HR_First, subject_number, c("lead_val", "SUB_NR"), sep="_")
df_HR_First$SUB_NR <- as.integer(df_HR_First$SUB_NR)
df_HR <- merge(df_HR, df_HR_First, by="SUB_NR", all=T)
df_HR$first_scan <- paste(df_HR$first_scan, "First", sep="-")
df_HR <- dplyr::filter(df_HR, !is.na(first_scan))
#Make Separate df for early and late
df_HR_Early <- subset(df_HR, time_point=='Early')
df_HR_Late <- subset(df_HR, time_point=='Late')
df_HR_Early$Run_time[df_HR_Early$Run_time==1] <- 6
# Clean names
df_HR <- df_HR %>% clean_names() %>% filter(first_scan!="NA-First")
df_HR$subject_number <- paste("sub", str_pad(df_HR$sub_nr, 3, pad = "0"), sep="_")
df_HR <- merge(df_HR, df_demographics, by="subject_number") %>% dplyr::rename(first_scan="first_scan.x")

2.2 Statistics

Now that we have cleaned up our physiology data, we can start with the actual statistics of interest. The goal here is to establish increased physiological arousal/stress in the stress session relative to the control. We first take a look at relevant saliva measures (i.e., cortisol and alpha amylase), and then the heart rate measures (using IBI and RMSSD).

Cortisol

We first check the data distribution, which looks like it needs a log transformation. We then run a model, with the addition of time, and controlling for scan order effects. We see in our main model that there is a significant stress*time effect on cortisol. We run a follow-up on the initial model to check the individual time effects. We see that that during the early stress phase, there is a significant increase increase in cortisol levels. We see that there is a significant difference in cortisol between the stress and control sessions, and at specific points in the follow up.

Since we want to use a single metric in the fMRI analysis for the response as a modulator, we check the area under the curve (AUC) between the two time points. We do this with a paired sample t-test since its the difference in a single measure we are testing here. We see there is a significant difference in the AUCi between the sessions. We also test an alternative method for baseline corrections (the Henckens response). To test if there is an increase in cortisol in the stress session, this value should be positive (stress-control/baseline). If stress-control is > 0, then the stress work. So we use a one-sample t-test against 0 to check if it works here.

Hist.

hist(df_saliva.scan$cortisol)

hist(log(df_saliva.scan$cortisol))

Main Effects

model.cort_scan <- glmer(cortisol ~ scan_type*time*first_scan + sex + (1+scan_type|subject_number), 
                         family=Gamma(link="log"),
                         glmerControl(calc.derivs =FALSE),
                         contrasts=list(scan_type="contr.treatment", sex="contr.treatment"),
                         data=df_saliva.scan)
asis_output(tab_model(model.cort_scan, show.stat=TRUE, show.se=TRUE, transform=NULL)$knitr)

	cortisol
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	1.43	0.06	1.30 – 1.55	22.28	<0.001
scan_type [Stress]	-0.03	0.07	-0.17 – 0.11	-0.45	0.654
time	-0.00	0.00	-0.00 – -0.00	-4.63	<0.001
first_scan1	0.05	0.06	-0.06 – 0.16	0.82	0.411
sex [Male]	-0.16	0.08	-0.32 – -0.01	-2.06	0.040
scan_type [Stress] * time	0.00	0.00	0.00 – 0.00	2.44	0.015
scan_type [Stress] * first_scan1	-0.03	0.07	-0.17 – 0.10	-0.50	0.618
time * first_scan1	0.00	0.00	-0.00 – 0.00	0.18	0.855
scan_type [Stress] * time * first_scan1	0.00	0.00	-0.00 – 0.00	0.07	0.941
Random Effects
σ2	0.18
τ00 subject_number	0.09
τ11 subject_number.scan_type1	0.05
ρ01 subject_number	0.06
ICC	0.33
N subject_number	76
Observations	713
Marginal R2 / Conditional R2	0.049 / 0.367

check_model(model.cort_scan)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Loading required namespace: qqplotr

Post-Hoc

df_saliva.scan$run <- as.factor(df_saliva.scan$time_point)
model.cort_scan_followup <- glmer(cortisol ~ scan_type*run + scan_type*first_scan + sex + (1+scan_type|subject_number), 
                                  contrasts=list(scan_type="contr.treatment"),
                         family=Gamma(link="log"), 
                         glmerControl(calc.derivs =FALSE),
                         data=df_saliva.scan)
emmeans(model.cort_scan_followup, list(pairwise ~ scan_type*run),by=c("run"), adjust="none")

$`emmeans of scan_type | run`
run = 0:
 scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control     1.33 0.0653 Inf     1.206      1.46
 Stress      1.21 0.0634 Inf     1.087      1.34

run = 14.2042253521127:
 scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control     1.36 0.0652 Inf     1.231      1.49
 Stress      1.19 0.0637 Inf     1.062      1.31

run = 42.7816901408451:
 scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control     1.25 0.0653 Inf     1.119      1.38
 Stress      1.54 0.0637 Inf     1.419      1.67

run = 87.0070422535211:
 scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control     1.09 0.0656 Inf     0.957      1.21
 Stress      1.25 0.0634 Inf     1.126      1.37

run = 159.724637681159:
 scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control     1.09 0.0656 Inf     0.966      1.22
 Stress      1.20 0.0634 Inf     1.071      1.32

Results are averaged over the levels of: first_scan, sex 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of scan_type | run`
run = 0:
 2                estimate     SE  df z.ratio p.value
 Control - Stress    0.123 0.0838 Inf   1.464  0.1431

run = 14.2042253521127:
 2                estimate     SE  df z.ratio p.value
 Control - Stress    0.172 0.0840 Inf   2.050  0.0404

run = 42.7816901408451:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   -0.296 0.0840 Inf  -3.529  0.0004

run = 87.0070422535211:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   -0.165 0.0840 Inf  -1.962  0.0498

run = 159.724637681159:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   -0.101 0.0839 Inf  -1.206  0.2280

Results are averaged over the levels of: first_scan, sex 
Results are given on the log (not the response) scale. 

Post-Hoc: Sex/Contraceptive

model.cort_scan_contr <- glmer(cortisol ~ scan_type*time*first_scan + contraceptive_use + (1+scan_type|subject_number), 
                         family=Gamma(link="log"),
                         glmerControl(calc.derivs =FALSE),
                         contrasts=list(scan_type="contr.treatment", contraceptive_use="contr.treatment"),
                         data=df_saliva.scan)
asis_output(tab_model(model.cort_scan_contr, show.stat=TRUE, show.se=TRUE, transform=NULL)$knitr)

	cortisol
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	1.25	0.07	1.10 – 1.40	16.79	<0.001
scan_type [Stress]	-0.01	0.07	-0.16 – 0.13	-0.20	0.844
time	-0.00	0.00	-0.00 – -0.00	-4.05	<0.001
first_scan1	0.07	0.06	-0.05 – 0.18	1.14	0.253
contraceptive_use [No]	0.15	0.10	-0.05 – 0.35	1.50	0.134
contraceptive_use [Yes]	0.26	0.10	0.07 – 0.45	2.69	0.007
scan_type [Stress] * time	0.00	0.00	0.00 – 0.00	2.29	0.022
scan_type [Stress] * first_scan1	-0.04	0.07	-0.19 – 0.10	-0.61	0.543
time * first_scan1	-0.00	0.00	-0.00 – 0.00	-0.33	0.739
scan_type [Stress] * time * first_scan1	0.00	0.00	-0.00 – 0.00	0.24	0.811
Random Effects
σ2	0.18
τ00 subject_number	0.09
τ11 subject_number.scan_type1	0.05
ρ01 subject_number	0.05
ICC	0.33
N subject_number	69
Observations	662
Marginal R2 / Conditional R2	0.067 / 0.376

check_model(model.cort_scan_contr)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

AUCi

df_saliva.wide_sub <- df_saliva.scan %>% dplyr::select(subject_number, age, sex, contraceptive_use, first_scan, AUCg_control:HPA_act) %>% unique() 
describe(df_saliva.wide$AUCi_difference<0)

df_saliva.wide$AUCi_difference < 0 
       n  missing distinct 
      64       23        2 
                      
Value      FALSE  TRUE
Frequency     50    14
Proportion 0.781 0.219

model.auc <- lm(AUCi_stress - AUCi_control ~ 1 + first_scan + sex, data=df_saliva.wide_sub)
summary(model.auc)

Call:
lm(formula = AUCi_stress - AUCi_control ~ 1 + first_scan + sex, 
    data = df_saliva.wide_sub)

Residuals:
    Min      1Q  Median      3Q     Max 
-503.67 -118.63   -4.56  114.89  662.41 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    95.20      25.03   3.804 0.000332 ***
first_scan1   -39.95      24.88  -1.606 0.113501    
sex1           69.29      25.08   2.763 0.007558 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 198.7 on 61 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.1506,    Adjusted R-squared:  0.1227 
F-statistic: 5.408 on 2 and 61 DF,  p-value: 0.006886

HPA Response

model.hpa <- lm(HPA_act ~ 1 + first_scan + sex, data=df_saliva.wide_sub)
summary(model.hpa)

Call:
lm(formula = HPA_act ~ 1 + first_scan + sex, data = df_saliva.wide_sub)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6825 -0.5174 -0.0295  0.4185  3.4922 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.18614    0.14896   1.250    0.216
first_scan1 -0.01271    0.14699  -0.086    0.931
sex1         0.11093    0.14914   0.744    0.460

Residual standard error: 1.228 on 67 degrees of freedom
  (9 observations deleted due to missingness)
Multiple R-squared:  0.008424,  Adjusted R-squared:  -0.02118 
F-statistic: 0.2846 on 2 and 67 DF,  p-value: 0.7532

Plot

# Make Summary
df_saliva.summ <- summarySEwithin(data=df_saliva.scan, idvar ="subject_number", withinvars =c("scan_type", "time_point"), measurevar = "cortisol", na.rm=TRUE)

Automatically converting the following non-factors to factors: time_point

df_saliva.summ$time_point <- as.character(df_saliva.summ$time_point) %>% as.numeric(df_saliva.summ$time_point)
#Construct the plot
plot.cort_time <- ggplot(df_saliva.summ, aes(y=cortisol, x=time_point,  colour=scan_type)) +
  # Add the scan phases
  geom_rect(aes(xmin=3, xmax=10, ymin=3, ymax=6.2, color=NULL), fill="#E68613", alpha=0.02) +
  geom_rect(aes(xmin=18, xmax=84, ymin=3, ymax=6.2, color=NULL), fill="lightgrey", alpha=0.02) +
  geom_rect(aes(xmin=91, xmax=157, ymin=3, ymax=6.2 ,  color=NULL),fill="lightgrey", alpha=0.02) +
  # Actial plot
  geom_line(size=1, stat = "identity") +
  geom_errorbar(aes(ymin=cortisol-se, ymax=cortisol+se),size=1, width=3, na.rm=T) +
  geom_point(shape=20, size=3) + 
  # Scale to make the graph start at 3
  scale_y_continuous(expand = c(0, 0), limits = c(3, NA)) +
  # Titles and labels
  labs(colour="Session", fill="Phase")+
  xlab("Time (minutes)") + ylab("Cortisol (nmol/L)\n") +
  ggtheme2 + scale_color_manual(values=stress_colors)
#Plot with time 
plot.cort_time; ggsave(filename="figures/Figure_3_LabStres_A.svg", device="svg", plot = plot.cort_time, dpi=300, bg="white") #

Amylase

In addition to the cortisol, we also assayed the levels of alpha amylase, which is used a proxy for adrenaline. We first take a look at the distribution of our variable to see if we need to do any transformations to it before modeling, and then proceed to do the models in the same way we did for the cortisol measures. We see that the distribution is a bit off, and that we need to apply a log transformation to get better fit in our models.

In the main effects of our models, we see that there is a significant effect of time, and that there is a significant interaction between scan order and the stress effects. We investigate this in follow-up models, and see that in people who had a stress session first, the heart rate was higher in the stress than control sessions, but not in those who had the control session first. We go on to plot this effect.

Hist.

#hist(df_saliva.scan$alpha_amylase)
hist(log(df_saliva.scan$alpha_amylase))

Main Effects

model.amylase_scan <- glmer(alpha_amylase ~ scan_type*time*first_scan + sex + (1+scan_type|subject_number), 
                            contrasts=list(scan_type="contr.treatment", sex="contr.treatment"),
                          family=Gamma(link="log"),
                          glmerControl(calc.derivs = FALSE),
                          data=df_saliva.scan)
asis_output(tab_model(model.amylase_scan, show.se=TRUE, show.stat=TRUE, transform=NULL)$knitr)

	alpha amylase
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	4.03	0.10	3.84 – 4.21	42.37	<0.001
scan_type [Stress]	-0.07	0.07	-0.22 – 0.07	-0.97	0.330
time	0.00	0.00	0.00 – 0.00	6.38	<0.001
first_scan1	-0.22	0.08	-0.38 – -0.07	-2.79	0.005
sex [Male]	0.29	0.13	0.03 – 0.55	2.19	0.029
scan_type [Stress] * time	0.00	0.00	-0.00 – 0.00	0.88	0.379
scan_type [Stress] * first_scan1	0.17	0.07	0.03 – 0.32	2.35	0.019
time * first_scan1	0.00	0.00	-0.00 – 0.00	0.14	0.889
scan_type [Stress] * time * first_scan1	0.00	0.00	-0.00 – 0.00	0.05	0.958
Random Effects
σ2	0.22
τ00 subject_number	0.29
τ11 subject_number.scan_type1	0.05
ρ01 subject_number	0.05
ICC	0.57
N subject_number	75
Observations	695
Marginal R2 / Conditional R2	0.129 / 0.623

check_model(model.amylase_scan)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

Post-Hoc

model.amylase_scan.fu <- glmer(alpha_amylase ~ scan_type*run + scan_type*first_scan + sex +(1+scan_type|subject_number), 
                               contrasts=list(scan_type="contr.treatment"),
                          family=Gamma(link="log"),
                          glmerControl(calc.derivs = FALSE),
                          data=df_saliva.scan)
emmeans(model.amylase_scan.fu, pairwise ~ first_scan*scan_type, adjust="none")

$emmeans
 first_scan scan_type emmean     SE  df asymp.LCL asymp.UCL
 Control    Control     4.12 0.1034 Inf      3.92      4.32
 Stress     Control     4.60 0.1010 Inf      4.40      4.79
 Control    Stress      4.24 0.1013 Inf      4.05      4.44
 Stress     Stress      4.34 0.0996 Inf      4.15      4.54

Results are averaged over the levels of: run, sex 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast                         estimate     SE  df z.ratio p.value
 Control Control - Stress Control  -0.4773 0.1433 Inf  -3.332  0.0009
 Control Control - Control Stress  -0.1255 0.0919 Inf  -1.365  0.1722
 Control Control - Stress Stress   -0.2233 0.1423 Inf  -1.569  0.1167
 Stress Control - Control Stress    0.3519 0.1417 Inf   2.483  0.0130
 Stress Control - Stress Stress     0.2540 0.0878 Inf   2.893  0.0038
 Control Stress - Stress Stress    -0.0978 0.1408 Inf  -0.695  0.4871

Results are averaged over the levels of: run, sex 
Results are given on the log (not the response) scale. 

#emmip(model.amylase_scan.fu, ~ scan_type*first_scan)

Plot

# Make Summary
df_saliva.summ_amy <-  summarySEwithin(data=df_saliva.scan, idvar ="subject_number", withinvars =c("scan_type", "time_point"), betweenvars = "first_scan", measurevar = "alpha_amylase", na.rm=TRUE)

Automatically converting the following non-factors to factors: first_scan, time_point

df_saliva.summ_amy$time_point <- as.character(df_saliva.summ_amy$time_point) %>% as.numeric(df_saliva.summ_amy$time_point)
#Construct the plot
plot.amylase_time <- ggplot(df_saliva.summ_amy , aes(y=alpha_amylase, x=time_point, colour=scan_type)) +
  geom_line(size=1,stat = "identity", na.rm = TRUE, show.legend = TRUE) +
  geom_errorbar(aes(ymin=alpha_amylase-se, ymax=alpha_amylase+se),size=1, width=3, na.rm=T) + 
  geom_point(shape=20, size=3) + 
  geom_rect(mapping=aes(xmin=3, xmax=12, ymin=50, ymax=150, fill='grey'), fill= 'lightgrey',color='lightgrey', alpha=0.05) +
  scale_y_continuous()+# scale_x_continuous(breaks=c(0, 6, 18, 50, 100))+
  ggtitle("Salivary Amylase During Scanning")+labs(colour="Session")+
  xlab("Time (minutes)") + ylab("Amylase (nmol/L)\n")+ 
  ggtheme2 + scale_color_manual(values=stress_colors)
#Plot with error bats and marking significant differences
plot.amylase_time + facet_grid(.~first_scan)

IBI

Now we can look at the changes in the IBI as a function of stress. We expect the IBI to be increased following the stress induction. Here, we analyze the effects of across the 10 scanner runs our participants had (including the resting state runs).In the models, we again include the scan order as a fixed effect in the same way we did for the cortisol analysis.

There is a significant main effect of time, a two-way interaction between stress and scan order, and three-way interaction between stress effects, time, and scan order. A follow-up shows that the stress effect is only present in those who have the stress session first,

Main Effects

# IBI Linear Model
model.IBI_scan <- lmer(ibi ~ session*run*first_scan + sex + (session|sub_nr),
                       contrasts=list(session="contr.treatment", sex="contr.treatment"),
                       data=df_HR)
asis_output(tab_model(model.IBI_scan, show.se=T, show.stat=T)$knitr)

	ibi
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	67.61	1.22	65.23 – 70.00	55.60	<0.001
session [Stress]	1.30	0.81	-0.29 – 2.89	1.60	0.110
run	-0.11	0.04	-0.18 – -0.04	-2.91	0.004
first_scan1	0.97	0.99	-0.96 – 2.90	0.98	0.325
sex [Male]	-4.70	1.84	-8.29 – -1.10	-2.56	0.011
session [Stress] * run	0.07	0.05	-0.04 – 0.17	1.24	0.213
session [Stress] * first_scan1	-2.50	0.81	-4.10 – -0.91	-3.08	0.002
run * first_scan1	-0.12	0.04	-0.20 – -0.05	-3.21	0.001
session [Stress] * run * first_scan1	0.20	0.05	0.10 – 0.31	3.83	<0.001
Random Effects
σ2	8.10
τ00 sub_nr	59.12
τ11 sub_nr.session1	9.99
ρ01 sub_nr	-0.02
ICC	0.88
N sub_nr	75
Observations	1410
Marginal R2 / Conditional R2	0.091 / 0.890

check_model(model.IBI_scan)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

Post-Hoc

df_HR$Run <- as.factor(df_HR$run)
model.ibi_scan.fu <- lmer(ibi ~ session*Run*first_scan + sex + (session|sub_nr), 
                           contrasts=list(session="contr.treatment", sex="contr.treatment"),
                         data=df_HR )
emmeans(model.ibi_scan.fu, list(pairwise ~ session*first_scan),by=c("first_scan"), adjust="none")

NOTE: Results may be misleading due to involvement in interactions

$`emmeans of session | first_scan`
first_scan = Control-First:
 session emmean   SE   df lower.CL upper.CL
 Control   64.9 1.38 71.6     62.2     67.7
 Stress    65.2 1.42 72.9     62.4     68.1

first_scan = Stress-First:
 session emmean   SE   df lower.CL upper.CL
 Control   64.4 1.38 72.6     61.6     67.1
 Stress    67.4 1.39 71.3     64.6     70.2

Results are averaged over the levels of: Run, sex 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of session | first_scan`
first_scan = Control-First:
 2                estimate   SE   df t.ratio p.value
 Control - Stress    -0.28 1.09 70.9  -0.258  0.7975

first_scan = Stress-First:
 2                estimate   SE   df t.ratio p.value
 Control - Stress    -3.04 1.06 70.6  -2.868  0.0054

Results are averaged over the levels of: Run, sex 
Degrees-of-freedom method: kenward-roger 

emmeans(model.ibi_scan.fu, list(pairwise ~ session*Run),by=c("first_scan", "Run"), adjust="none")

$`emmeans of session | first_scan, Run`
first_scan = Control-First, Run = 1:
 session emmean   SE   df lower.CL upper.CL
 Control   67.7 1.44 84.4     64.8     70.5
 Stress    66.0 1.48 85.6     63.1     68.9

first_scan = Control-First, Run = 2:
 session emmean   SE   df lower.CL upper.CL
 Control   64.6 1.43 83.5     61.7     67.4
 Stress    63.9 1.47 85.1     60.9     66.8

first_scan = Control-First, Run = 3:
 session emmean   SE   df lower.CL upper.CL
 Control   65.8 1.43 83.5     62.9     68.6
 Stress    65.9 1.47 85.1     63.0     68.9

first_scan = Control-First, Run = 4:
 session emmean   SE   df lower.CL upper.CL
 Control   64.9 1.43 83.5     62.1     67.8
 Stress    64.9 1.48 85.6     62.0     67.8

first_scan = Control-First, Run = 5:
 session emmean   SE   df lower.CL upper.CL
 Control   66.7 1.43 83.5     63.9     69.6
 Stress    66.7 1.47 85.1     63.7     69.6

first_scan = Control-First, Run = 6:
 session emmean   SE   df lower.CL upper.CL
 Control   63.2 1.44 84.5     60.4     66.1
 Stress    64.3 1.48 86.1     61.3     67.2

first_scan = Control-First, Run = 7:
 session emmean   SE   df lower.CL upper.CL
 Control   62.1 1.44 85.5     59.2     65.0
 Stress    63.8 1.48 86.6     60.8     66.7

first_scan = Control-First, Run = 8:
 session emmean   SE   df lower.CL upper.CL
 Control   64.7 1.43 84.0     61.8     67.5
 Stress    65.2 1.48 86.1     62.2     68.1

first_scan = Control-First, Run = 9:
 session emmean   SE   df lower.CL upper.CL
 Control   64.1 1.44 84.5     61.2     67.0
 Stress    64.4 1.48 85.6     61.5     67.3

first_scan = Control-First, Run = 10:
 session emmean   SE   df lower.CL upper.CL
 Control   65.7 1.44 84.9     62.8     68.5
 Stress    67.3 1.48 85.6     64.3     70.2

first_scan = Stress-First, Run = 1:
 session emmean   SE   df lower.CL upper.CL
 Control   65.8 1.44 85.1     62.9     68.6
 Stress    69.8 1.45 83.5     66.9     72.6

first_scan = Stress-First, Run = 2:
 session emmean   SE   df lower.CL upper.CL
 Control   63.8 1.44 85.1     60.9     66.6
 Stress    66.2 1.45 84.4     63.3     69.1

first_scan = Stress-First, Run = 3:
 session emmean   SE   df lower.CL upper.CL
 Control   64.6 1.44 85.1     61.8     67.5
 Stress    68.2 1.45 84.5     65.3     71.0

first_scan = Stress-First, Run = 4:
 session emmean   SE   df lower.CL upper.CL
 Control   63.5 1.43 84.7     60.6     66.3
 Stress    67.0 1.45 83.6     64.1     69.9

first_scan = Stress-First, Run = 5:
 session emmean   SE   df lower.CL upper.CL
 Control   65.5 1.44 85.6     62.7     68.4
 Stress    69.4 1.45 84.4     66.5     72.3

first_scan = Stress-First, Run = 6:
 session emmean   SE   df lower.CL upper.CL
 Control   63.1 1.44 85.6     60.2     65.9
 Stress    66.4 1.44 82.7     63.5     69.2

first_scan = Stress-First, Run = 7:
 session emmean   SE   df lower.CL upper.CL
 Control   62.3 1.44 85.1     59.5     65.2
 Stress    64.7 1.44 82.7     61.9     67.6

first_scan = Stress-First, Run = 8:
 session emmean   SE   df lower.CL upper.CL
 Control   64.7 1.44 85.6     61.9     67.6
 Stress    67.4 1.44 83.1     64.5     70.3

first_scan = Stress-First, Run = 9:
 session emmean   SE   df lower.CL upper.CL
 Control   64.1 1.44 86.6     61.2     67.0
 Stress    66.1 1.44 83.1     63.2     69.0

first_scan = Stress-First, Run = 10:
 session emmean   SE   df lower.CL upper.CL
 Control   66.4 1.44 86.1     63.5     69.2
 Stress    69.0 1.44 82.7     66.1     71.9

Results are averaged over the levels of: sex 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of session | first_scan, Run`
first_scan = Control-First, Run = 1:
 3                estimate   SE  df t.ratio p.value
 Control - Stress  1.69983 1.23 116   1.380  0.1702

first_scan = Control-First, Run = 2:
 3                estimate   SE  df t.ratio p.value
 Control - Stress  0.70438 1.22 113   0.575  0.5663

first_scan = Control-First, Run = 3:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -0.14584 1.22 113  -0.119  0.9054

first_scan = Control-First, Run = 4:
 3                estimate   SE  df t.ratio p.value
 Control - Stress  0.00501 1.23 114   0.004  0.9967

first_scan = Control-First, Run = 5:
 3                estimate   SE  df t.ratio p.value
 Control - Stress  0.04284 1.22 113   0.035  0.9722

first_scan = Control-First, Run = 6:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -1.01956 1.23 116  -0.826  0.4104

first_scan = Control-First, Run = 7:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -1.69517 1.24 120  -1.365  0.1749

first_scan = Control-First, Run = 8:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -0.48492 1.23 116  -0.394  0.6946

first_scan = Control-First, Run = 9:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -0.30529 1.23 116  -0.248  0.8047

first_scan = Control-First, Run = 10:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -1.60524 1.23 117  -1.301  0.1960

first_scan = Stress-First, Run = 1:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -3.98740 1.20 116  -3.320  0.0012

first_scan = Stress-First, Run = 2:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -2.43760 1.21 117  -2.021  0.0455

first_scan = Stress-First, Run = 3:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -3.55847 1.21 117  -2.951  0.0038

first_scan = Stress-First, Run = 4:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -3.53374 1.20 115  -2.945  0.0039

first_scan = Stress-First, Run = 5:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -3.84504 1.21 118  -3.183  0.0019

first_scan = Stress-First, Run = 6:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -3.30250 1.20 115  -2.754  0.0068

first_scan = Stress-First, Run = 7:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -2.40852 1.20 114  -2.013  0.0465

first_scan = Stress-First, Run = 8:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -2.69218 1.20 116  -2.241  0.0269

first_scan = Stress-First, Run = 9:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -2.00561 1.21 117  -1.663  0.0990

first_scan = Stress-First, Run = 10:
 3                estimate   SE  df t.ratio p.value
 Control - Stress -2.64410 1.20 116  -2.200  0.0298

Results are averaged over the levels of: sex 
Degrees-of-freedom method: kenward-roger 

Plot

# Make Summary
df_HR.summ <- summarySEwithin(data=df_HR, idvar ="sub_nr", withinvars =c("session", "run"), betweenvars = "first_scan", measurevar = "ibi", na.rm=TRUE)

Automatically converting the following non-factors to factors: first_scan, run
NaNs producedNaNs produced

df_HR.summ <- na.omit(df_HR.summ)
df_HR.summ$run <- as.numeric(as.character(df_HR.summ$run))
df_HR.summ$Phase <- ifelse(df_HR.summ$run<6, "Early", "Late")
#Construct the plot
plot.ibi_time <- ggplot(df_HR.summ, aes(y=ibi, x=run, color=session)) +
  # Make the background session squares
  geom_rect(aes(xmin=0, xmax=0.6, ymin=62, ymax=72, color=NULL), fill="#E68613", alpha=0.02) +
  geom_rect(aes(xmin=0.8, xmax=5.2, ymin=62, ymax=72, color=NULL), fill="lightgrey", alpha=0.02) +
  geom_rect(aes(xmin=5.8, xmax=10.2, ymin=62, ymax=72 ,  color=NULL),fill="lightgrey", alpha=0.02) +
  # Actual Plot
  geom_line(stat = "identity", show.legend = TRUE) +
  geom_errorbar(aes(ymin=ibi-se, ymax=ibi+se), width=0.25, na.rm=T) +
  geom_point(shape=20, size=2, line=2) + 
  # Legen titles and such
  labs(colour="Session", fill="Session")+
  xlab("") + ylab("Heart Rate (BPM)")+  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10)) + 
  guides(color=guide_legend(order=2)) +
  ggtheme2 + scale_color_manual(values=stress_colors )  #+ scale_fill_manual(values=c("lightgrey", "lightgrey", ))

Ignoring unknown parameters: line

#plot.ibi_time
plot.ibi_time + facet_grid(first_scan~.)

RMSSD

Now we can take a look at the RMSSD values. First we check the distribution to aid with model decisions and then run the models. There is a significant time by session effect here, with a significant effect of scan order on the session effect. We see that in people who had a stress session first, there is an overall lower RMSSD in the stress session.

Hist.

hist(log(df_HR$rtmssd))

Main Effects

model.rmssd_scan <- glmer(rtmssd ~ session*run*first_scan + sex + (session|sub_nr),
                           contrasts=list(session="contr.treatment", sex="contr.treatment"),
                         data=df_HR,
                         family=Gamma(link="log"),
                         control=glmerControl(calc.derivs=FALSE))
asis_output(tab_model(model.rmssd_scan, show.stat=TRUE, show.se=TRUE, transform=NULL, digits =3 )$knitr)

	rtmssd
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	3.847	0.046	3.758 – 3.937	84.438	<0.001
session [Stress]	0.019	0.042	-0.064 – 0.102	0.458	0.647
run	0.015	0.003	0.008 – 0.022	4.386	<0.001
first_scan1	-0.020	0.039	-0.097 – 0.056	-0.518	0.604
sex [Male]	0.049	0.060	-0.070 – 0.167	0.803	0.422
session [Stress] * run	-0.016	0.005	-0.025 – -0.006	-3.246	0.001
session [Stress] * first_scan1	0.118	0.042	0.036 – 0.201	2.800	0.005
run * first_scan1	0.006	0.003	-0.001 – 0.013	1.791	0.073
session [Stress] * run * first_scan1	-0.012	0.005	-0.022 – -0.003	-2.518	0.012
Random Effects
σ2	0.07
τ00 sub_nr	0.06
τ11 sub_nr.session1	0.02
ρ01 sub_nr	0.05
ICC	0.47
N sub_nr	75
Observations	1410
Marginal R2 / Conditional R2	0.039 / 0.492

check_model(model.rmssd_scan)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

Post-Hoc

df_HR$Run <- as.factor(df_HR$run)
model.rmssd_scan.fu <- glmer(rtmssd ~ session*Run*first_scan + sex + (session|sub_nr), 
                              contrasts=list(session="contr.treatment"),
                         data=df_HR,
                         family=Gamma(link="log"),
                         control=glmerControl(calc.derivs=FALSE, optimizer ="bobyqa"))

maxfun < 10 * length(par)^2 is not recommended.

emmeans(model.rmssd_scan.fu, list(pairwise ~ session*first_scan),by=c("first_scan"), adjust="none")

NOTE: Results may be misleading due to involvement in interactions

$`emmeans of session | first_scan`
first_scan = Control-First:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.97 0.0482 Inf      3.87      4.06
 Stress    3.95 0.0474 Inf      3.86      4.05

first_scan = Stress-First:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.94 0.0483 Inf      3.85      4.04
 Stress    3.82 0.0463 Inf      3.73      3.91

Results are averaged over the levels of: Run, sex 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of session | first_scan`
first_scan = Control-First:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.0162 0.0463 Inf   0.350  0.7265

first_scan = Stress-First:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.1201 0.0453 Inf   2.650  0.0081

Results are averaged over the levels of: Run, sex 
Results are given on the log (not the response) scale. 

emmeans(model.rmssd_scan.fu, list(pairwise ~ session*Run),by=c("Run"), adjust="none")

NOTE: Results may be misleading due to involvement in interactions

$`emmeans of session | Run`
Run = 1:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.85 0.0448 Inf      3.76      3.94
 Stress    3.88 0.0442 Inf      3.79      3.96

Run = 2:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.91 0.0445 Inf      3.83      4.00
 Stress    3.92 0.0444 Inf      3.83      4.00

Run = 3:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.91 0.0445 Inf      3.82      4.00
 Stress    3.85 0.0444 Inf      3.76      3.93

Run = 4:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.93 0.0444 Inf      3.84      4.02
 Stress    3.89 0.0443 Inf      3.80      3.98

Run = 5:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.90 0.0447 Inf      3.81      3.99
 Stress    3.85 0.0444 Inf      3.76      3.94

Run = 6:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   4.07 0.0450 Inf      3.98      4.15
 Stress    3.91 0.0442 Inf      3.82      4.00

Run = 7:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   4.02 0.0452 Inf      3.93      4.11
 Stress    3.95 0.0443 Inf      3.86      4.04

Run = 8:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.98 0.0449 Inf      3.89      4.07
 Stress    3.92 0.0443 Inf      3.83      4.01

Run = 9:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   4.01 0.0453 Inf      3.92      4.10
 Stress    3.89 0.0441 Inf      3.81      3.98

Run = 10:
 session emmean     SE  df asymp.LCL asymp.UCL
 Control   3.97 0.0453 Inf      3.88      4.06
 Stress    3.82 0.0440 Inf      3.73      3.90

Results are averaged over the levels of: first_scan, sex 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of session | Run`
Run = 1:
 2                estimate     SE  df z.ratio p.value
 Control - Stress -0.02615 0.0521 Inf  -0.502  0.6156

Run = 2:
 2                estimate     SE  df z.ratio p.value
 Control - Stress -0.00162 0.0520 Inf  -0.031  0.9751

Run = 3:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.06477 0.0520 Inf   1.247  0.2125

Run = 4:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.03769 0.0517 Inf   0.729  0.4663

Run = 5:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.04939 0.0521 Inf   0.948  0.3429

Run = 6:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.15717 0.0521 Inf   3.015  0.0026

Run = 7:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.06973 0.0525 Inf   1.329  0.1838

Run = 8:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.06161 0.0521 Inf   1.182  0.2372

Run = 9:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.11743 0.0524 Inf   2.242  0.0250

Run = 10:
 2                estimate     SE  df z.ratio p.value
 Control - Stress  0.15117 0.0523 Inf   2.892  0.0038

Results are averaged over the levels of: first_scan, sex 
Results are given on the log (not the response) scale. 

Plot

# Make Summary
df_HR.summ <- summarySEwithin(data=df_HR, idvar ="sub_nr", withinvars =c("session", "run"),betweenvars ="first_scan", measurevar = "rtmssd", na.rm=TRUE)

Automatically converting the following non-factors to factors: first_scan, run
NaNs producedNaNs produced

df_HR.summ <- na.omit(df_HR.summ)
df_HR.summ$run <- as.numeric(as.character(df_HR.summ$run))
#Construct the plot
plot.rtmssd_time <- ggplot(df_HR.summ, aes(y=rtmssd, x=run, colour=session)) +
  # Make the background session squares
  geom_rect(aes(xmin=0, xmax=0.6, ymin=47, ymax=72, color=NULL), fill="#E68613", alpha=0.02) +
  geom_rect(aes(xmin=0.8, xmax=5.2, ymin=47, ymax=72, color=NULL), fill="lightgrey", alpha=0.02) +
  geom_rect(aes(xmin=5.8, xmax=10.2, ymin=47, ymax=72 ,  color=NULL),fill="lightgrey", alpha=0.02) +
  # Add labels
  # geom_label(aes(x=0.4, y=75, label="SECPT"), color="black") +
  # geom_label(aes(x=3, y=75, label="Early"), color="black") +
  # geom_label(aes(x=8, y=75, label="Late"), color="black") +
  # Actual Plot
  geom_line(stat = "identity", show.legend = TRUE) +
  geom_errorbar(aes(ymin=rtmssd-se, ymax=rtmssd+se), width=0.25, na.rm=T) +
  geom_point(shape=20, size=2, line=2) + 
  labs(colour="Session", fill="Session")+
  xlab("Scanner Run") + ylab("RMSSD (ms)") +  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9,10)) +scale_y_continuous(n.breaks = 5) +
  ggtheme2 + scale_color_manual(values=stress_colors)

Ignoring unknown parameters: line

plot.rtmssd_time + facet_grid(first_scan~.) 

Full Plot

plot.physio_full <- ggarrange(plot.cort_time, 
                              ggarrange(NULL, (plot.ibi_time + theme(axis.text.x=element_text(size=6)) + facet_grid(first_scan~.)), 
                                       NULL, (plot.rtmssd_time + theme(axis.text.x=element_text(size=6))+ facet_grid(first_scan~.)),
                                       ncol=2, nrow=2, legend="none", align = "hv", labels = c("B", NA, "C", NA), widths=c(0.1,1)), 
          common.legend = T, legend="bottom", widths=c(1.5, 1), labels=c("A"))

# I separated this chunk because ggarrange was producing a random blank plot
plot.physio_full
ggsave("figures/Figure_3_LabStress.svg", dpi=300, device="svg", bg="white", width =8)

Saving 8 x 4.5 in image

3 fMRI: Large-scale networks under stress

We have now validated our stress exposure paradigms in both the lab and real life. We were also able to derive a single measure from each of the stress exposure paradigms we have. That is, we have the AUCi as our indicator of lab stress reactivity, and our residual based stress score for the real-life data. We can incorporate these two measures into our fMRI analysis. We do this in two ways: GLMs in the FSL software tool box, and models here in R for our ROI based analysis.

Before we start this analysis though, we do some signal exploration below. We describe briefly in the Introduction what tasks were used, but not the entirety of the fMRI design. In summary, our participants underwent a stress and control fMRI session. These sessions consisted of two halves, each with 3 runs. The three runs in each half are modeled at the first level in FSL to end up with four fMRI sessions per participant. The first three runs consist of the early phase, while the last three consist of the late phase. This results in:

• Control Early
• Control Late
• Stress Early
• Stress Late

This design gives us insight into the temporal dynamics of stress reactivity and recovery. We first explore the mean task related activity across all four sessions, and then examine the stress effects (stress-control). These effects are modeled at the subject level in FSL. For both the mean task and the stress related differences, we have two contrasts. The baseline contrast models the mean task activity relative to the implicit baseline, while the task related contrasts models the task related activity during the specific task, relative to the activity during the other tasks. This latter contrast is the main one we are interested in, as it specifically reflects the network resource allocations.

# Load data
df_fmri.signal <-read.csv("/project/3013068.02/stats/fMRI/Signal/fMRI_signal_event_stress_zstat.txt", sep='\t', header=T) %>%
  dplyr::group_by(sub_nr, cope) %>% distinct(Mask,Run, .keep_all=TRUE) %>% ungroup %>% dplyr::mutate(Network=Mask) %>% 
  dplyr::mutate(Contrast=str_replace(cope, "cope1", "2Back"),                                                                     
                Contrast=str_replace(Contrast, "cope2", "NoBack"),
                Contrast=str_replace(Contrast, "cope3", "2Back>NoBack"), 
                Contrast=str_replace(Contrast, "cope4", "Oddball"),
                Contrast=str_replace(Contrast, "cope5", "Noddball"),
                Contrast=str_replace(Contrast, "cope6", "Odd>Nodd"),
                Contrast=str_replace(Contrast, "cope7", "Mem-Rem"),
                Contrast=str_replace(Contrast, "cope8", "Mem-For"),
                Contrast=str_replace(Contrast, "cope9", "Rem>For"),
                Contrast=str_replace(Contrast, "2Back0", "ECN-Con"),
                Contrast=str_replace(Contrast, "2Back1", "SN-Con"),
                Contrast=str_replace(Contrast, "2Back2", "DMN-Con"),
                Contrast=str_replace(Contrast, "2Back3", "2Back+NoBack")) %>%
    dplyr::mutate(Contrast= factor(Contrast, levels=c("2Back", "NoBack", "2Back>NoBack", "ECN-Con", "2Back+NoBack",
                                                      "Oddball", "Noddball", "Odd>Nodd","SN-Con", 
                                                      "Mem-Rem", "Mem-For", "Rem>For", "DMN-Con")))
# Filter for plots and average tests
df_fmri.signal_mean <- df_fmri.signal %>% filter(Run=="Mean")
df_fmri.signal_event <- df_fmri.signal %>% filter(Run!="Mean")

We will also try to link the laboratory stress measure, to the real-life measure. We then try to see how either of these could be related to our large scale network balance. In order to do so, we need to perform some data cleaning, where we add all the data of interest together. We need to combine the fMRI, the cortisol, and the stress-based residual together into one data frame. We also save a data frame with just three columns: subject number, cortisol, and the residualized score to use in FSL directly.

# Merge the dataframes to include the stress reactivity measures
df_combo <- merge(df_stressrx, df_saliva.wide_sub, by.x="id", by.y="subject_number") 
df_combo.fmri <- df_fmri.signal_event
df_combo.fmri$id <- paste("sub", (str_pad(df_combo.fmri$sub_nr, 3, pad = "0")), sep="_") 
df_combo <- merge(df_combo, df_combo.fmri, by="id")
# Rescale and center data points for modeling
## Cortisol
df_combo$HPA_act[is.na(df_combo$HPA_act)] <- mean(df_combo$HPA_act, na.rm=T)
df_combo$AUCi_difference[is.na(df_combo$AUCi_difference)] <- mean(df_combo$AUCi_difference, na.rm=T)
## Rescale
df_combo <- df_combo %>% subset(! is.na(Signal)) %>% dplyr::mutate(HPA_act_z=scale(HPA_act, center=TRUE), stress_rx_z=scale(resBIN, center=TRUE), signal_z=scale(Signal, center=TRUE), AUCi_diff_z=scale(AUCi_difference, center=TRUE))
df_combo <- as.data.frame(df_combo)
# Contrast settings
df_combo$sex <- as.factor(df_combo$sex)
contrasts(df_combo$sex) <- contr.treatment(c("Female", "Male"))
df_combo$contraceptive_use <- as.factor(df_combo$contraceptive_use)
contrasts(df_combo$contraceptive_use) <- contr.treatment(c("Male", "No", "Yes"))
## Save as dataframe for fMRI whole brain analysis
df_combo.tosave <- df_combo %>% dplyr::select(sub_nr, stress_rx_z, HPA_act_z, AUCi_diff_z) %>% unique()
write.csv2(df_combo.tosave, "data/fMRI_regressor.txt", dec="." )

attempt to set 'dec' ignored

3.1 Mean Task Effects

Here we check for the main task effects for our 2Back (ECN), Oddball (SN), and Retrieval (DMN) contrasts. To do this, we plot the mean signal from each the tasks across networks. What we would ideally like to see is the recruitment of the corresponding network-task combinations. While we do see this in the ECN and SN contrasts, we do not really get this effect for the DMN.

ECN: 2Back

# Filter and generate summary
df_fmri.2back <- df_combo %>% subset(cope=="cope1" | cope=="cope2" | cope=="cope3")
df_fmri.2back <- summarySEwithin(data=df_fmri.2back, idvar="sub_nr", measurevar="Signal", withinvars = c("Network", "cope", "Contrast"))
asis_output(tab_model(lmer(Signal ~ 1 + sex+ (1|id), data=(df_combo%>% subset(cope=="cope3" & Network=="ECN"))))$knitr)

	Signal
Predictors	Estimates	CI	p
(Intercept)	0.57	0.41 – 0.74	<0.001
sex [Male]	-0.00	-0.27 – 0.26	0.974
Random Effects
σ2	0.46
τ00 id	0.19
ICC	0.29
N id	70
Observations	280
Marginal R2 / Conditional R2	0.000 / 0.292

# Make plot
plot_fmri.mean_ecn <- ggplot(df_fmri.2back , aes(y=Signal, x=Network, color=Network, fill=Network, na.rm = TRUE)) +
  geom_bar(stat="identity", alpha=0.5, position="dodge2" ) +
   geom_errorbar(aes(ymin=Signal-se, ymax=Signal+se, color=Network),size=1, width=0.5, na.rm=T, position=position_dodge(.9)) +
  geom_hline(yintercept=0, color="grey") +
  ggtitle("Parameter Estimates of Task-Network Engagement") + ylab("P.E.") + xlab(" ") + 
  scale_color_manual(values=fmri_colors) + scale_fill_manual(values=fmri_colors) + 
  ggtheme2 + theme(axis.text.y=element_text(size=12), plot.title=element_text(hjust = 0), legend.position="right")
plot_fmri.mean_ecn  + facet_wrap(Contrast~., ncol=3)

SN: Oddball

# Subset + Model
df_fmri.odd <- df_combo %>% subset(cope=="cope4" | cope=="cope5" | cope=="cope6")
df_fmri.odd <- summarySEwithin(data=df_fmri.odd, idvar="sub_nr", measurevar="Signal", withinvars = c("Network", "cope", "Contrast"))
asis_output(tab_model(lmer(Signal ~ 1 + sex + (1|id), data=(df_combo%>% subset(cope=="cope6" & Network=="SN"))))$knitr)

	Signal
Predictors	Estimates	CI	p
(Intercept)	1.09	0.94 – 1.25	<0.001
sex [Male]	0.04	-0.20 – 0.29	0.742
Random Effects
σ2	0.57
τ00 id	0.12
ICC	0.17
N id	70
Observations	280
Marginal R2 / Conditional R2	0.001 / 0.169

# Plot
plot_fmri.mean_sn <- ggplot(df_fmri.odd, aes(y=Signal, x=Network, color=Network, fill=Network, na.rm = TRUE)) +
  geom_bar(stat="identity", alpha=0.5, position="dodge2" ) +
   geom_errorbar(aes(ymin=Signal-se, ymax=Signal+se, color=Network),size=1, width=0.5, na.rm=T, position=position_dodge(.9)) +
  geom_hline(yintercept=0, color="grey") +
  ggtitle("Parameter Estimates of Task-Network Engagement") + ylab("P.E.") + xlab(" ")+ 
  scale_color_manual(values=fmri_colors) + scale_fill_manual(values=fmri_colors) + 
  ggtheme2 + theme(axis.text.y=element_text(size=12), plot.title=element_text(hjust = 0), legend.position="right")
plot_fmri.mean_sn   + facet_wrap(Contrast~.)

DMN: Memory

# Subset + Model
df_fmri.mem <- df_combo %>% subset(cope=="cope7" | cope=="cope8" | cope=="cope9") 
df_fmri.mem <- summarySEwithin(data=df_fmri.mem, idvar="sub_nr", measurevar="Signal", withinvars = c("Network", "cope", "Contrast"))
asis_output(tab_model(lmer(Signal ~ 1 + sex + (1|id), data=(df_combo%>% subset(cope=="cope9" & Network=="DMN"))))$knitr)

	Signal
Predictors	Estimates	CI	p
(Intercept)	-0.40	-0.52 – -0.28	<0.001
sex [Male]	0.13	-0.06 – 0.33	0.176
Random Effects
σ2	0.53
τ00 id	0.03
ICC	0.06
N id	70
Observations	280
Marginal R2 / Conditional R2	0.008 / 0.065

# Plot
plot_fmri.mean_dmn <- ggplot(df_fmri.mem, aes(y=Signal, x=Network, color=Network, fill=Network, na.rm = TRUE)) +
  geom_bar(stat="identity", alpha=0.5, position="dodge2" ) +
   geom_errorbar(aes(ymin=Signal-se, ymax=Signal+se, color=Network),size=1, width=0.5, na.rm=T, position=position_dodge(.9)) +
  geom_hline(yintercept=0, color="grey") +
  ggtitle("Parameter Estimates of Task-Network Engagement") + ylab("P.E.") + xlab(" ") + 
  scale_color_manual(values=fmri_colors) + scale_fill_manual(values=fmri_colors) + 
  ggtheme2 + theme(axis.text.y=element_text(size=12), plot.title=element_text(hjust = 0), legend.position="right")
plot_fmri.mean_dmn   + facet_wrap(Contrast~.)

ROI-Task Pairs

df_fmri.mean <- df_combo %>% subset((cope=="cope3" & Network=="ECN") | (cope=="cope6" & Network=="SN") | (cope=="cope9" & Network=="DMN"))
df_fmri.mean_sum <- summarySEwithin(data=df_fmri.mean, idvar="sub_nr", measurevar="Signal", withinvars = c("Network", "cope", "Contrast"))
df_fmri.mean_sum$Contrast <- factor(df_fmri.mean_sum$Contrast)
plot.fmri_mean <- ggplot(df_fmri.mean_sum, aes(y=Signal, x=Network, color=Contrast, fill=Contrast, na.rm = TRUE)) +
  geom_bar(stat="identity", alpha=0.5, position="dodge2", width=0.7 ) +
   geom_errorbar(aes(ymin=Signal-se, ymax=Signal+se),size=1, width=0.4, na.rm=T, position=position_dodge(.9)) +
  geom_hline(yintercept=0, color="grey") +
  ggtitle("Parameter Estimates of Task-Network Engagement") + ylab("Parameter Estiamtes") + xlab("\nROI") + 
  scale_color_manual(values=c("#0CB702", "#F8766D", "#00BFC4")) + scale_fill_manual(values=c("#0CB702", "#F8766D", "#00BFC4" )) + 
  ggtheme2 + 
  theme(plot.title=element_text(hjust = 0), legend.position="right", 
        axis.text.y=element_text(size=12), 
        axis.text.x=element_text(size=12, colour="black"))
plot.fmri_mean

3.2 Stress Related Activity

We next take a look at the difference between stress and control sessions in the MRI. We first code the stress and control session as well as the early and late phases to make some factors out of them. We also filter only the specific contrasts we want to investigate. That is, the 2Back > non-target, the oddball > standard, and the remembered > forgotten contrasts. Once we have that, we construct a first model to look at main effects, and then break down these effects per task/network. We also can take a look at the difference in neural stress reactivity, and how that can further inform us regarding daily life stress reactivity. Finally we take a look at the concept of network balance, as laid out by Krause et. al. 2021.

df_fmri.stress <- df_combo %>% separate(Run, into=c("Session","Run"), sep="-" ) %>% subset((cope=="cope3" & Network=="ECN") | (cope=="cope6" & Network=="SN") | (cope=="cope9" & Network=="DMN")) %>% dplyr::mutate(Network=relevel(Network, "SN", "ECN", "DMN"))

Main Effects

In our main model, we see that there is a significant session by ROI and session by run interaction effect. The mixed models show us a comparison relative to baseline conditions (Control, Early) and in network 1 which we coded as SN. We see that there is a significant network effect, session*run interaction, and session by network interaction. We also see there is a significant 4-way interaction effect between session, run, network, and the resilience score in the ECN relative to SN. To simplify the interpretation, we also present the models in a simple ANOVA output for easier interpretation for the model with only the SN and ECN. The establishment of these main effects also allows us to run the next follow-up tests.

Model: Full

# Main model without covaritates
model.fmri_main <- lmer(Signal ~ Session*Run*Network + first_scan + sex + (1+Session+Run+Network|id), 
                         contrasts=list(Session="contr.treatment"),
                        data=df_fmri.stress, 
                        control=lmerControl(optimizer="bobyqa", calc.derivs = F))
# Residal Score Model
model.fmri_main_res <- lmer(Signal ~ Session*Run*Network*resBIN + first_scan +sex + (1+Session+Run+Network|id), 
                            contrasts=list(Session="contr.treatment"),
                        data=df_fmri.stress, 
                        control=lmerControl(optimizer="bobyqa", calc.derivs = F))
# AUCi Model
model.fmri_main_auci <- lmer(Signal ~ Session*Run*Network*AUCi_difference + first_scan + sex +(1+Session+Run+Network|id), 
                             contrasts=list(Session="contr.treatment"),
                        data=df_fmri.stress, 
                        control=lmerControl(optimizer="bobyqa", calc.derivs = F))
# Full Model
model.fmri_main_full <- lmer(Signal ~ Session*Run*Network*resBIN + Session*Run*Network*AUCi_difference + first_scan + sex + (1+Session+Run+Network|id), 
                             contrasts=list(Session="contr.treatment"),
                        data=df_fmri.stress, 
                        control=lmerControl(optimizer="bobyqa", calc.derivs = F))
# HTML Table with all
asis_output(tab_model(model.fmri_main,  model.fmri_main_res, model.fmri_main_auci,  model.fmri_main_full,
                             show.stat=TRUE, show.se=TRUE, digits = 3, 
                             dv.labels =c("Signal", "Signal~Resilience", "Signal~AUCi",  "Signal ~Full"))$knitr)

	Signal					Signal~Resilience					Signal~AUCi					Signal ~Full
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.426	0.054	0.320 – 0.531	7.882	<0.001	0.431	0.054	0.326 – 0.537	7.986	<0.001	0.361	0.063	0.237 – 0.486	5.708	<0.001	0.371	0.064	0.245 – 0.497	5.771	<0.001
Session [Stress]	-0.035	0.054	-0.141 – 0.071	-0.642	0.521	-0.035	0.054	-0.142 – 0.072	-0.641	0.521	0.000	0.063	-0.123 – 0.123	0.003	0.998	0.000	0.063	-0.124 – 0.124	0.003	0.997
Run1	-0.013	0.034	-0.081 – 0.054	-0.379	0.705	-0.013	0.034	-0.081 – 0.054	-0.387	0.699	-0.026	0.040	-0.105 – 0.053	-0.649	0.516	-0.026	0.040	-0.105 – 0.052	-0.650	0.516
Network1	0.589	0.058	0.474 – 0.703	10.077	<0.001	0.586	0.059	0.471 – 0.701	9.994	<0.001	0.593	0.068	0.459 – 0.727	8.679	<0.001	0.593	0.069	0.459 – 0.727	8.656	<0.001
Network2	-0.706	0.056	-0.816 – -0.596	-12.590	<0.001	-0.701	0.056	-0.810 – -0.592	-12.590	<0.001	-0.727	0.065	-0.855 – -0.598	-11.099	<0.001	-0.727	0.065	-0.855 – -0.600	-11.185	<0.001
first_scan1	-0.058	0.035	-0.126 – 0.011	-1.659	0.097	-0.061	0.035	-0.129 – 0.008	-1.723	0.085	-0.046	0.035	-0.115 – 0.023	-1.303	0.192	-0.049	0.036	-0.119 – 0.021	-1.374	0.170
sex [Male]	0.084	0.072	-0.056 – 0.225	1.174	0.240	0.063	0.074	-0.082 – 0.208	0.847	0.397	0.115	0.074	-0.031 – 0.261	1.545	0.122	0.095	0.078	-0.057 – 0.248	1.224	0.221
Session [Stress] * Run1	0.114	0.049	0.019 – 0.209	2.344	0.019	0.115	0.048	0.020 – 0.210	2.366	0.018	0.137	0.057	0.027 – 0.248	2.431	0.015	0.137	0.056	0.027 – 0.248	2.437	0.015
Session [Stress] * Network1	0.156	0.069	0.021 – 0.290	2.266	0.023	0.159	0.069	0.025 – 0.294	2.327	0.020	0.140	0.080	-0.016 – 0.297	1.755	0.079	0.140	0.080	-0.016 – 0.296	1.755	0.079
Session [Stress] * Network2	-0.175	0.069	-0.310 – -0.040	-2.549	0.011	-0.179	0.069	-0.313 – -0.045	-2.612	0.009	-0.112	0.080	-0.269 – 0.045	-1.397	0.162	-0.111	0.080	-0.268 – 0.045	-1.396	0.163
Run1 * Network1	-0.077	0.049	-0.172 – 0.018	-1.585	0.113	-0.081	0.048	-0.176 – 0.014	-1.678	0.093	-0.075	0.057	-0.185 – 0.036	-1.319	0.187	-0.074	0.056	-0.185 – 0.036	-1.314	0.189
Run1 * Network2	0.017	0.049	-0.078 – 0.112	0.346	0.729	0.019	0.048	-0.076 – 0.114	0.391	0.696	0.042	0.057	-0.069 – 0.152	0.735	0.462	0.041	0.056	-0.069 – 0.152	0.733	0.463
Session [Stress] * Run1 * Network1	0.103	0.069	-0.032 – 0.238	1.501	0.133	0.108	0.069	-0.026 – 0.243	1.582	0.114	0.085	0.080	-0.071 – 0.242	1.069	0.285	0.085	0.080	-0.071 – 0.241	1.064	0.287
Session [Stress] * Run1 * Network2	-0.030	0.069	-0.164 – 0.105	-0.431	0.667	-0.030	0.069	-0.164 – 0.104	-0.436	0.663	-0.065	0.080	-0.222 – 0.092	-0.811	0.417	-0.065	0.080	-0.221 – 0.091	-0.814	0.416
resBIN						0.085	0.066	-0.044 – 0.214	1.291	0.197						0.069	0.065	-0.059 – 0.197	1.057	0.290
Session [Stress] * resBIN						0.008	0.076	-0.141 – 0.157	0.105	0.916						0.014	0.076	-0.135 – 0.163	0.184	0.854
Run1 * resBIN						0.009	0.048	-0.085 – 0.103	0.184	0.854						0.007	0.048	-0.088 – 0.101	0.138	0.890
Network1 * resBIN						0.093	0.082	-0.067 – 0.254	1.143	0.253						0.095	0.082	-0.067 – 0.256	1.148	0.251
Network2 * resBIN						-0.166	0.078	-0.318 – -0.014	-2.135	0.033						-0.170	0.078	-0.323 – -0.017	-2.175	0.030
(Session [Stress] * Run1) * resBIN						-0.030	0.068	-0.162 – 0.102	-0.442	0.658						-0.026	0.068	-0.159 – 0.107	-0.383	0.702
(Session [Stress] Network1) resBIN						-0.140	0.095	-0.327 – 0.047	-1.467	0.142						-0.143	0.096	-0.331 – 0.045	-1.495	0.135
(Session [Stress] Network2) resBIN						0.144	0.095	-0.043 – 0.332	1.513	0.130						0.156	0.096	-0.032 – 0.344	1.627	0.104
(Run1 * Network1) * resBIN						0.157	0.068	0.025 – 0.290	2.329	0.020						0.158	0.068	0.026 – 0.291	2.337	0.019
(Run1 * Network2) * resBIN						-0.078	0.068	-0.210 – 0.054	-1.154	0.248						-0.074	0.068	-0.207 – 0.059	-1.093	0.274
(Session [Stress] * Run1 * Network1) * resBIN						-0.196	0.095	-0.383 – -0.009	-2.054	0.040						-0.200	0.096	-0.388 – -0.012	-2.087	0.037
(Session [Stress] * Run1 * Network2) * resBIN						0.011	0.095	-0.176 – 0.198	0.113	0.910						0.005	0.096	-0.183 – 0.193	0.051	0.960
AUCi_difference											0.000	0.000	-0.000 – 0.001	1.888	0.059	0.000	0.000	-0.000 – 0.001	1.725	0.084
Session [Stress] * AUCi_difference											-0.000	0.000	-0.001 – 0.000	-1.090	0.276	-0.000	0.000	-0.001 – 0.000	-1.093	0.274
Run1 * AUCi_difference											0.000	0.000	-0.000 – 0.000	0.635	0.525	0.000	0.000	-0.000 – 0.000	0.625	0.532
Network1 * AUCi_difference											-0.000	0.000	-0.001 – 0.001	-0.113	0.910	-0.000	0.000	-0.001 – 0.001	-0.196	0.845
Network2 * AUCi_difference											0.000	0.000	-0.000 – 0.001	0.630	0.529	0.000	0.000	-0.000 – 0.001	0.789	0.430
(Session [Stress] * Run1) * AUCi_difference											-0.000	0.000	-0.001 – 0.000	-0.822	0.411	-0.000	0.000	-0.001 – 0.000	-0.795	0.427
(Session [Stress] Network1) AUCi_difference											0.000	0.000	-0.001 – 0.001	0.371	0.711	0.000	0.000	-0.001 – 0.001	0.479	0.632
(Session [Stress] Network2) AUCi_difference											-0.001	0.000	-0.001 – 0.000	-1.549	0.121	-0.001	0.000	-0.001 – 0.000	-1.667	0.095
(Run1 * Network1) * AUCi_difference											-0.000	0.000	-0.001 – 0.000	-0.081	0.936	-0.000	0.000	-0.001 – 0.000	-0.250	0.803
(Run1 * Network2) * AUCi_difference											-0.000	0.000	-0.001 – 0.000	-0.859	0.390	-0.000	0.000	-0.001 – 0.000	-0.781	0.435
(Session [Stress] * Run1 * Network1) * AUCi_difference											0.000	0.000	-0.001 – 0.001	0.430	0.667	0.000	0.000	-0.000 – 0.001	0.581	0.561
(Session [Stress] * Run1 * Network2) * AUCi_difference											0.000	0.000	-0.000 – 0.001	0.866	0.386	0.000	0.000	-0.000 – 0.001	0.863	0.388
Random Effects
σ2	0.49					0.49					0.50					0.49
τ00	0.06 id					0.06 id					0.06 id					0.06 id
τ11	0.01 id.Session1					0.01 id.Session1					0.01 id.Session1					0.01 id.Session1
	0.00 id.Run1					0.00 id.Run1					0.00 id.Run1					0.00 id.Run1
	0.07 id.Network1					0.08 id.Network1					0.08 id.Network1					0.08 id.Network1
	0.06 id.Network2					0.05 id.Network2					0.06 id.Network2					0.05 id.Network2
ρ01	-0.11					-0.11					-0.20					-0.19
	0.10					0.14					0.07					0.10
	0.00					-0.01					-0.00					-0.02
	-0.71					-0.68					-0.70					-0.67
N	70 id					70 id					70 id					70 id
Observations	840					840					840					840
Marginal R2 / Conditional R2	0.435 / NA					0.444 / NA					0.438 / NA					0.447 / NA

Model: Full as ANOVA

asis_output(kbl(anova(model.fmri_main), digits=3, caption ="**Classic ANOVA Output of Main Effect**", align ="l"))

Classic ANOVA Output of Main Effect

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
Session	0.204	0.204	1	69.678	0.413	0.523
Run	1.586	1.586	1	431.460	3.206	0.074
Network	168.901	84.450	2	77.450	170.672	0.000
first_scan	1.362	1.362	1	82.180	2.753	0.101
sex	0.682	0.682	1	82.180	1.379	0.244
Session:Run	2.719	2.719	1	619.446	5.495	0.019
Session:Network	3.863	1.931	2	619.445	3.903	0.021
Run:Network	0.336	0.168	2	619.446	0.339	0.713
Session:Run:Network	1.182	0.591	2	619.446	1.194	0.304

Plot

# Calculate the difference between stress and control 
df_fmri.stress_wide <- df_fmri.stress %>% pivot_wider(., id_cols = c("sub_nr","age","sex", "contraceptive_use", "first_scan", "Network","Contrast"), names_from =c("Session","Run" ), values_from=c("Signal", "resBIN", "AUCi_difference")) %>% 
  mutate(`SE-CE`=Signal_Stress_Early-Signal_Control_Early, `SL-CL`=Signal_Stress_Late-Signal_Control_Late) %>%
  pivot_longer(., cols=c( "SE-CE", "SL-CL"), names_to="Run",values_to="Stress-Control" ) %>% 
  dplyr::mutate(Run=str_replace(Run, "SE-CE", "Early"), Run=str_replace(Run, "SL-CL", "Late"))
df_fmri.stress$sex <- as.factor(df_fmri.stress$sex)
contrasts(df_fmri.stress$sex) <- contr.treatment(c("Male", "Female"))
# Make summary stats
df_fmri.stress_summ <- summarySEwithin(data=df_fmri.stress_wide, idvar="sub_nr", measurevar="Stress-Control", withinvars = c("Run", "Network", "Contrast"))
df_fmri.stress_summ$Net_Run <-paste(df_fmri.stress_summ$Network, df_fmri.stress_summ$Run, sep="\n")
df_fmri.stress_summ$Net_Run <- fct_relevel(df_fmri.stress_summ$Net_Run, "DMN-Early", "DMN-Late", "ECN-Early", "ECN-Late", "SN-Early")
# Plot
fmri_colors2 <- c("#00BFC4","#00B8E7", "#7CAE00", "#0CB702", "#F8766D","#E68613" )
plot.fmri_stress <- ggplot(df_fmri.stress_summ, aes(y=`Stress-Control`, x=Net_Run, color=Net_Run, fill=Net_Run, na.rm = TRUE)) +
   geom_hline(yintercept=0, color="grey") + 
   geom_bar(stat="identity", alpha=0.5, position=position_dodge2() ) +
   geom_errorbar(aes(ymin=`Stress-Control`-se, ymax=`Stress-Control`+se),size=1, na.rm=T,width=0.4, position=position_dodge2(padding = 0.75)) +
  ylab(sprintf("\U03B2Stress - \U03B2 Control (a.u.)")) + xlab("Networks") +
   scale_color_manual(values=fmri_colors2) + scale_fill_manual(values=fmri_colors2) + 
   ggtheme2 + theme(plot.title=element_text(hjust=0, size=16), legend.position="none",
                    axis.text.y=element_text(size=10), axis.title.y=element_text(size=12),
                    axis.text.x=element_text(size=10), axis.ticks.x=element_line(size=1), axis.title.x =element_text(size=12)) 
  #guides( fill="none",color=guide_legend(colour=fmri_colors, override.aes=list(linetype=c(1, 1, 1), shape=c(16,16,16), fill=c("#00BFC4","#0CB702","#F8766D") ))) 
plot.fmri_stress; ggsave("figures/Figure_4_A.svg", device="svg", bg="white")

Network Effects

We established some interesting stress effects above. Now we want to perform follow-up tests for each network individually. Note that these results to some extent (for example, session and run effects) are just repetitions of the main effects we established above. Interestingly, we see that there is a significant interaction between stress, time, and the residual resilience score in the salience network.

ECN

# Filter df
df_fmri.ecn <- df_fmri.stress %>% subset(Contrast=="2Back>NoBack" & Network=="ECN" )
# Run model with AUCI
model.ecn_auc <- lmer(Signal ~ Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                      contrasts=list(Session="contr.treatment"),
                      data=df_fmri.ecn, 
                      control=lmerControl(optimizer=c("bobyqa"), calc.derivs=F))
# Run model with resilience
model.ecn_res <- lmer(Signal ~ Session*Run*resBIN + first_scan + sex + (1+Session+Run|sub_nr),
                      contrasts=list(Session="contr.treatment"),
                      data=df_fmri.ecn, 
                      control=lmerControl(optimizer=c("bobyqa"), calc.derivs=F))
# Full model
model.ecn_full <- lmer(Signal ~ Session*Run*resBIN + Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                       contrasts=list(Session="contr.treatment"),
                      data=df_fmri.ecn,
                      control=lmerControl(optimizer=c("bobyqa"), calc.derivs=F))
# Tab and print results
tab.ecn_react <- tab_model(model.ecn_auc, model.ecn_res, model.ecn_full,
                           dv.labels=c("AUCi", "Real Life", "Full"), 
                           show.se = T, show.stat = T, emph.p = T,
                           df.method="kenward" )
knitr::asis_output(tab.ecn_react$knitr)

	AUCi					Real Life					Full
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.56	0.09	0.38 – 0.74	6.08	<0.001	0.60	0.09	0.42 – 0.77	6.56	<0.001	0.58	0.09	0.40 – 0.76	6.26	<0.001
Session [Stress]	-0.02	0.08	-0.18 – 0.15	-0.18	0.857	-0.02	0.09	-0.18 – 0.15	-0.18	0.857	-0.02	0.09	-0.18 – 0.15	-0.18	0.858
Run1	0.05	0.06	-0.06 – 0.16	0.85	0.398	0.05	0.06	-0.06 – 0.16	0.88	0.380	0.05	0.06	-0.06 – 0.16	0.88	0.377
AUCi_diff_z	0.07	0.08	-0.08 – 0.22	0.90	0.370						0.05	0.08	-0.10 – 0.20	0.64	0.523
first_scan1	0.08	0.07	-0.05 – 0.21	1.14	0.255	0.05	0.07	-0.08 – 0.17	0.72	0.471	0.06	0.07	-0.07 – 0.19	0.88	0.379
sex [Female]	0.05	0.14	-0.23 – 0.32	0.34	0.732	-0.06	0.14	-0.33 – 0.21	-0.43	0.667	-0.02	0.15	-0.30 – 0.27	-0.12	0.901
Session [Stress] * Run1	0.04	0.08	-0.11 – 0.19	0.51	0.608	0.04	0.08	-0.12 – 0.19	0.46	0.648	0.04	0.08	-0.12 – 0.19	0.46	0.647
Session [Stress] * AUCi_diff_z	0.02	0.08	-0.14 – 0.19	0.26	0.794						0.02	0.09	-0.15 – 0.19	0.26	0.796
Run1 * AUCi_diff_z	0.07	0.06	-0.04 – 0.18	1.22	0.224						0.07	0.06	-0.04 – 0.18	1.29	0.198
(Session [Stress] * Run1) * AUCi_diff_z	-0.13	0.08	-0.28 – 0.02	-1.64	0.101						-0.14	0.08	-0.29 – 0.02	-1.77	0.077
resBIN						0.16	0.11	-0.05 – 0.37	1.50	0.134	0.15	0.11	-0.07 – 0.36	1.35	0.178
Session [Stress] * resBIN						0.00	0.12	-0.23 – 0.24	0.03	0.976	0.00	0.12	-0.23 – 0.24	0.01	0.991
Run1 * resBIN						-0.07	0.08	-0.22 – 0.08	-0.90	0.366	-0.08	0.08	-0.23 – 0.07	-1.00	0.318
(Session [Stress] * Run1) * resBIN						0.16	0.11	-0.06 – 0.37	1.41	0.158	0.17	0.11	-0.04 – 0.38	1.55	0.120
Random Effects
σ2	0.43					0.44					0.43
τ00	0.20 sub_nr					0.19 sub_nr					0.19 sub_nr
τ11	0.02 sub_nr.Session1					0.02 sub_nr.Session1					0.02 sub_nr.Session1
	0.00 sub_nr.Run1					0.00 sub_nr.Run1					0.00 sub_nr.Run1
ρ01	-0.34					-0.31					-0.29
	0.23					0.44					0.05
ICC	0.31										0.31
N	70 sub_nr					70 sub_nr					70 sub_nr
Observations	280					280					280
Marginal R2 / Conditional R2	0.030 / 0.331					0.053 / NA					0.050 / 0.349

SN

# Filter only SN
df_fmri.sn <- df_fmri.stress %>% subset(Contrast=="Odd>Nodd" & Network=="SN" )
# AUC model
model.sn_auc <- lmer(Signal ~ Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                     contrasts=list(Session="contr.treatment"),
                     data=df_fmri.sn)
# Resilience Model
model.sn_res <- lmer(Signal ~ Session*Run*resBIN + first_scan + sex + (1+Session+Run|sub_nr),
                     contrasts=list(Session="contr.treatment"),
                     data=df_fmri.sn)
# Full
model.sn_full <- lmer(Signal ~ Session*Run*resBIN + Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                      contrasts=list(Session="contr.treatment"),
                     data=df_fmri.sn)
# HTML Tab
tab.sn_react <- tab_model(model.sn_auc, model.sn_res, model.sn_full,
                          dv.labels=c("AUCi", "Real Life", "Full"), show.se = T, show.stat = T)
knitr::asis_output(tab.sn_react$knitr)

	AUCi					Real Life					Full
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	1.02	0.09	0.83 – 1.20	10.82	<0.001	1.05	0.09	0.86 – 1.23	11.31	<0.001	1.03	0.09	0.85 – 1.22	10.95	<0.001
Session [Stress]	0.12	0.10	-0.08 – 0.32	1.17	0.241	0.12	0.10	-0.08 – 0.33	1.21	0.225	0.12	0.10	-0.08 – 0.33	1.20	0.229
Run1	-0.09	0.06	-0.20 – 0.02	-1.54	0.124	-0.09	0.06	-0.21 – 0.02	-1.64	0.100	-0.09	0.06	-0.21 – 0.02	-1.63	0.103
AUCi_diff_z	0.08	0.08	-0.08 – 0.24	1.03	0.304						0.07	0.08	-0.09 – 0.23	0.82	0.413
first_scan1	0.00	0.06	-0.12 – 0.12	0.03	0.977	-0.02	0.06	-0.14 – 0.10	-0.30	0.766	-0.01	0.06	-0.13 – 0.11	-0.13	0.898
sex [Female]	0.08	0.13	-0.18 – 0.34	0.59	0.553	-0.01	0.13	-0.26 – 0.25	-0.06	0.951	0.03	0.14	-0.24 – 0.30	0.22	0.827
Session [Stress] * Run1	0.22	0.08	0.07 – 0.37	2.84	0.005	0.22	0.07	0.08 – 0.37	3.02	0.003	0.22	0.07	0.08 – 0.37	2.99	0.003
Session [Stress] * AUCi_diff_z	-0.03	0.10	-0.23 – 0.17	-0.32	0.747						-0.03	0.10	-0.23 – 0.18	-0.26	0.797
Run1 * AUCi_diff_z	0.02	0.06	-0.10 – 0.13	0.31	0.758						0.01	0.06	-0.10 – 0.12	0.16	0.871
(Session [Stress] * Run1) * AUCi_diff_z	-0.01	0.08	-0.16 – 0.14	-0.14	0.892						0.00	0.07	-0.14 – 0.15	0.02	0.986
resBIN						0.18	0.11	-0.04 – 0.40	1.64	0.101	0.17	0.11	-0.05 – 0.39	1.48	0.138
Session [Stress] * resBIN						-0.13	0.14	-0.41 – 0.15	-0.92	0.356	-0.13	0.14	-0.41 – 0.15	-0.90	0.370
Run1 * resBIN						0.17	0.08	0.01 – 0.32	2.07	0.038	0.17	0.08	0.01 – 0.32	2.04	0.041
(Session [Stress] * Run1) * resBIN						-0.23	0.10	-0.43 – -0.02	-2.19	0.028	-0.23	0.10	-0.43 – -0.02	-2.17	0.030
Random Effects
σ2	0.41					0.38					0.39
τ00	0.16 sub_nr					0.16 sub_nr					0.16 sub_nr
τ11	0.08 sub_nr.Session1					0.09 sub_nr.Session1					0.09 sub_nr.Session1
	0.04 sub_nr.Run1					0.04 sub_nr.Run1					0.04 sub_nr.Run1
ρ01	-0.03					-0.06					-0.07
	0.03					0.03					0.01
ICC	0.32					0.35					0.35
N	70 sub_nr					70 sub_nr					70 sub_nr
Observations	280					280					280
Marginal R2 / Conditional R2	0.033 / 0.346					0.053 / 0.384					0.057 / 0.383

emtrends(model.sn_full, identity ~ Session | Run, var="resBIN")

$emtrends
Run = Early:
 Session resBIN.trend    SE   df lower.CL upper.CL
 Control      0.33364 0.135 92.5   0.0664    0.601
 Stress      -0.02180 0.148 88.3  -0.3157    0.272

Run = Late:
 Session resBIN.trend    SE   df lower.CL upper.CL
 Control      0.00336 0.144 89.3  -0.2837    0.290
 Stress       0.10005 0.137 91.5  -0.1730    0.373

Results are averaged over the levels of: first_scan, sex 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
Run = Early:
 contrast estimate    SE   df t.ratio p.value
 Control   0.33364 0.135 92.5   2.479  0.0150
 Stress   -0.02180 0.148 88.3  -0.147  0.8832

Run = Late:
 contrast estimate    SE   df t.ratio p.value
 Control   0.00336 0.144 89.3   0.023  0.9815
 Stress    0.10005 0.137 91.5   0.728  0.4686

Results are averaged over the levels of: first_scan, sex 
Degrees-of-freedom method: kenward-roger 

DMN

# Filter DMN tasks
df_fmri.dmn <- df_fmri.stress %>% subset(Contrast=="Rem>For" & Network=="DMN" )
# AUC Model
model.dmn_auc <- lmer(Signal ~ Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                      contrasts=list(Session="contr.treatment"),
                      data=df_fmri.dmn)

Model failed to converge with max|grad| = 0.00369878 (tol = 0.002, component 1)

# Resilience Model
model.dmn_res <- lmer(Signal ~ Session*Run*resBIN + first_scan + sex + (1+Session+Run|sub_nr),
                      contrasts=list(Session="contr.treatment"),
                      data=df_fmri.dmn)
# Full model
model.dmn_full <- lmer(Signal ~ Session*Run*resBIN + Session*Run*AUCi_diff_z + first_scan + sex + (1+Session+Run|sub_nr),
                       contrasts=list(Session="contr.treatment"),
                       control=lmerControl(optimizer="bobyqa", calc.derivs = F),
                      data=df_fmri.dmn)
# HTML Output
tab.dmn_react <- tab_model(model.dmn_auc, model.dmn_res, model.dmn_full, 
                           dv.labels=c("AUCi", "Real Life", "Full"),show.se=T, show.est=T, show.stat = T)
knitr::asis_output(tab.dmn_react$knitr)

	AUCi					Real Life					Full
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	-0.31	0.08	-0.46 – -0.16	-4.00	<0.001	-0.30	0.08	-0.45 – -0.15	-3.84	<0.001	-0.31	0.08	-0.46 – -0.15	-3.94	<0.001
Session [Stress]	-0.21	0.09	-0.39 – -0.03	-2.34	0.019	-0.21	0.09	-0.39 – -0.04	-2.35	0.019	-0.21	0.09	-0.39 – -0.04	-2.41	0.016
Run1	0.00	0.06	-0.11 – 0.12	0.06	0.950	0.01	0.06	-0.11 – 0.12	0.10	0.924	0.01	0.06	-0.11 – 0.12	0.09	0.925
AUCi_diff_z	0.12	0.07	-0.01 – 0.25	1.75	0.080						0.13	0.07	-0.01 – 0.26	1.82	0.069
first_scan1	-0.10	0.05	-0.20 – -0.01	-2.18	0.029	-0.11	0.05	-0.20 – -0.02	-2.32	0.020	-0.10	0.05	-0.20 – -0.01	-2.12	0.034
sex [Female]	0.15	0.10	-0.04 – 0.35	1.53	0.126	0.13	0.10	-0.06 – 0.33	1.35	0.177	0.16	0.11	-0.05 – 0.37	1.49	0.137
Session [Stress] * Run1	0.08	0.08	-0.08 – 0.25	1.02	0.306	0.08	0.08	-0.08 – 0.25	1.03	0.303	0.08	0.08	-0.08 – 0.25	1.02	0.306
Session [Stress] * AUCi_diff_z	-0.16	0.09	-0.34 – 0.01	-1.85	0.065						-0.17	0.09	-0.35 – 0.00	-1.95	0.051
Run1 * AUCi_diff_z	-0.02	0.06	-0.14 – 0.10	-0.34	0.737						-0.02	0.06	-0.13 – 0.10	-0.28	0.783
(Session [Stress] * Run1) * AUCi_diff_z	0.02	0.08	-0.14 – 0.18	0.24	0.811						0.02	0.08	-0.14 – 0.18	0.25	0.803
resBIN						-0.08	0.10	-0.27 – 0.10	-0.88	0.379	-0.10	0.10	-0.29 – 0.09	-1.07	0.284
Session [Stress] * resBIN						0.15	0.13	-0.10 – 0.40	1.20	0.229	0.17	0.12	-0.07 – 0.41	1.36	0.172
Run1 * resBIN						-0.07	0.08	-0.23 – 0.09	-0.84	0.400	-0.07	0.08	-0.23 – 0.10	-0.81	0.416
(Session [Stress] * Run1) * resBIN						-0.02	0.11	-0.24 – 0.21	-0.17	0.868	-0.02	0.12	-0.25 – 0.21	-0.18	0.855
Random Effects
σ2	0.47					0.47					0.48
τ00	0.04 sub_nr					0.04 sub_nr					0.04 sub_nr
τ11	0.02 sub_nr.Session1					0.03 sub_nr.Session1					0.02 sub_nr.Session1
	0.01 sub_nr.Run1					0.01 sub_nr.Run1					0.01 sub_nr.Run1
ρ01	0.28					0.34					0.28
	0.49					0.54					0.55
ICC	0.09					0.08
N	70 sub_nr					70 sub_nr					70 sub_nr
Observations	280					280					280
Marginal R2 / Conditional R2	0.071 / 0.152					0.068 / 0.146					0.088 / NA

SN Reactivity

We saw the main effects above, now we want to look as specific follow-ups of the stress resilience scores and the changes in MRI stress activity in our networks over time. We will only do this for the SN, since that is the only one of our main tests where we see an effect of stress on the network, and an interaction effect with the residualized score. To this end, we subtract Stress-Control from the MRI sessions to get a single measure of within subject reactivity in the lab in the early MRI run, and recover in the late MRI run.

Additionally, in these models we slightly change the contrasts to see specific differences between the early and late phases. We therefore code the effect of Run as a contrast treatment here.

df_fmri.stress_wide <- df_fmri.stress %>%  pivot_wider(id_cols=c("id","age","sex","contraceptive_use", "Run", "Network", "first_scan"), names_from="Session",  values_from=c("Signal", "resBIN", "AUCi_difference")) %>%  dplyr::mutate(Signal=Signal_Stress-Signal_Control) %>% dplyr::rename(resBIN=resBIN_Control)

Main Model

df_fmri.stress_sn <- df_fmri.stress_wide %>% dplyr::filter(Network=="SN")
model.sn_run <- lmer(Signal ~ Run*resBIN + (1|id),  data=df_fmri.stress_sn, contrasts=list(Run="contr.treatment"))
asis_output(tab_model(model.sn_run, show.se=T, show.stat=TRUE)$knitr)

	Signal
Predictors	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.34	0.13	0.10 – 0.59	2.73	0.006
Run [Late]	-0.44	0.15	-0.72 – -0.15	-2.99	0.003
resBIN	-0.36	0.17	-0.71 – -0.02	-2.08	0.037
Run [Late] * resBIN	0.47	0.20	0.07 – 0.86	2.29	0.022
Random Effects
σ2	0.75
τ00 id	0.36
ICC	0.32
N id	70
Observations	142
Marginal R2 / Conditional R2	0.068 / 0.369

Post-Hoc

df_fmri.stress_sn_early <- df_fmri.stress_sn %>% filter(Run=="Early")
model.sn_run_early <- lm(Signal ~ resBIN, data=df_fmri.stress_sn_early)
df_fmri.stress_sn_late <- df_fmri.stress_sn %>% filter(Run=="Late")
model.sn_run_late <- lm(Signal ~  resBIN, data=df_fmri.stress_sn_late)
asis_output(tab_model(model.sn_run_early, model.sn_run_late, dv.labels = c("Early", "Late"), show.stat=TRUE, show.se=T)$knitr)

	Early					Late
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.34	0.12	0.10 – 0.59	2.77	0.007	-0.09	0.13	-0.34 – 0.16	-0.73	0.468
resBIN	-0.36	0.17	-0.71 – -0.02	-2.10	0.039	0.10	0.18	-0.25 – 0.45	0.58	0.563
Observations	71					71
R2 / R2 adjusted	0.060 / 0.047					0.005 / -0.010

emmip(model.sn_run, Run ~ resBIN, cov.reduce=range) + ggtheme2

Plot: Early SN-Resilience

# Plot the SN-Stress relaiton to Resilience
df_fmri.sn_early_wide <- df_fmri.sn %>% dplyr::filter(Run=="Early") %>% pivot_wider(id_cols=c("id", "age", "sex", "contraceptive_use", "first_scan"), names_from="Session",  values_from=c("resBIN", "Signal")) %>%  dplyr::mutate(Signal=Signal_Stress-Signal_Control) %>% dplyr::rename(resBIN=resBIN_Control) %>% dplyr::mutate(resBIN=resBIN*-1)
# Plot 1: Resilience by reactivity
plot.fmri_sndiff = ggplot(df_fmri.sn_early_wide, aes(x=resBIN, y=Signal, color="Stress-Control")) + 
  geom_smooth(method="lm",formula=y~x, se=F, color="#F8766D" ) +
  geom_point(alpha=0.7, color="#F8766D" ) + 
  xlab(sprintf("SN Early (\U03B2 Stress - \U03B2 Control)\n")) + ylab("Real Life Stressor Reactivity (a.u.)") +
  ggtheme2 + theme(  axis.text.y=element_text(size=10), axis.title.y=element_text(size=12),
                    axis.text.x=element_text(size=10), axis.ticks.x=element_line(size=1), axis.title.x =element_text(size=12))
plot.fmri_sndiff
ggsave("figures/Figure_4_B_SN.svg", device="svg", bg="white")

Saving 7.29 x 4.5 in image

Post-Hoc: Recovery

# Subset the data to get early-late difference
df_fmri.sn_wide_recover <- df_fmri.sn %>%  pivot_wider(id_cols=c("id", "age", "sex", "contraceptive_use", "first_scan"), names_from=c("Session","Run"),  values_from=c("resBIN", "Signal")) %>% 
  dplyr::mutate(Signal_Late=Signal_Stress_Late-Signal_Control_Late, Signal_Early=Signal_Stress_Early-Signal_Control_Early, Signal=Signal_Late-Signal_Early) %>%
  dplyr::rename(resBIN=resBIN_Stress_Early) 
# Remove outliers
signal_sd <- sd(df_fmri.sn_wide_recover$Signal*3)
signal_sd_neg <- -signal_sd
df_fmri.sn_wide_recover <- df_fmri.sn_wide_recover %>% mutate(Signal=replace(Signal, Signal<signal_sd_neg | Signal>signal_sd, NA))
# Model
asis_output(tab_model(lm(Signal ~ resBIN, data=df_fmri.sn_wide_recover), title="Change from Early to Late")$knitr)

Change from Early to Late

	Signal
Predictors	Estimates	CI	p
(Intercept)	-0.38	-0.66 – -0.10	0.008
resBIN	0.39	0.00 – 0.78	0.049
Observations	70
R2 / R2 adjusted	0.056 / 0.042

`

Full MRI Stress Plot

ggarrange(plot.fmri_stress, NULL, plot.fmri_sndiff, widths=c(1.2,0.01, 1), ncol = 3, labels = c("A", "B"), align = "hv")
ggsave("figures/Figure_4_fMRIStress.svg", device="svg")

Saving 7.29 x 4.5 in image

4 fMRI: Tasks

We have three tasks during the fMRI session which we need to analyze for our study as mentioned in the Introduction section. These tasks include a 2-Back, Oddball, and Associative Retrieval. Each requires a slightly different approach so we do some data cleaning in chunks of code before we start.

4.1 Data Cleaning

We first import the data from individual logs. To this end, we have a small function that reads a file and assigns the subject ID. We compile the trial-by-trial logs for all subjects into a single data frame which we then subset for each task separately. We then combine this data frame with the last data frame we used in the previous section for the balance.

# First get a list of all the log files
logs_beh <-(Sys.glob("/project/3013068.02/data/*/logs/mri/*mri*"))
# Make a function to read the log files + add session + subject info
read_add <- function(file, seperator, headers){
    df <- read.csv2(file=file, sep=seperator, header=headers)
    id <- file; id <- sub(".*sub_", "", id); id <- sub("/logs/.*", "", id)
    df$id <- id
    # Select session indicatior
    ses <- file; ses <- sub(".*logs/mri/Sub", "", ses); ses <- sub(".txt.*", "", ses)
    df$ses <- ses
    df }
# Run function
df_task <- rbindlist(lapply(logs_beh, read_add, seperator="\t", headers=FALSE))
# Coloumn headers
names(df_task) <- c("task", "stimulus", "valence", "trial_onset", "trial_offset", "rt", "button_pressed", "button_correct",
                         "block_onset","block_offset","experiment_start", "experiment_end", "sub_id", "session_id")
# Drop the headers added when importing, and retype
df_task <- df_task %>% subset(task != "Trial_Type" & task != "Trial_Type " ) %>% separate(session_id, c("sub", "mri", "session","run"))  %>% retype()
# Make a factor for task types and sessions and runs
df_task$sub_id <- as.factor(df_task$sub_id)
df_task$sub_nr <- paste("sub", (str_pad(df_task$sub_id, 3, pad = "0")), sep="_") 
# Label the tasks
df_task$Task <- factor(df_task$task, level=c(1,2,3), labels=c("oddball", "nback", "memory"))
# Label the session and runs
df_task$Session <- factor(df_task$session, level=c("s1","s2"), labels=c("Control", "Stress"))
df_task$Run <- factor(df_task$run, level=c("r1","r2","r3","r4","r5","r6"), labels=c("r1","r2","r3","r4","r5","r6"))
df_task$run <- as.numeric(df_task$Run)
df_task$Phase <- "Late"
df_task$Phase[df_task$Run =="r1" | df_task$Run =="r2" | df_task$Run =="r3"] <- "Early"
df_task$run <- as.numeric(df_task$Run)
# Subset the fmri dataframe
df_fmri.balance_4tasks <- df_fmri.stress  %>%  pivot_wider(id_cols=c(id, age, sex, contraceptive_use, first_scan, Session, Run), names_from=c(Network), values_from=c(Signal,resBIN,AUCi_difference) ) %>% dplyr::select(1:11,15) 
names(df_fmri.balance_4tasks ) <- c("sub_nr","age","sex", "contraceptive_use", "Scan_Order",  "Session", "Run", "ECN", "SN", "DMN", "resBIN", "AUCi_diff")
# Make the scores
df_fmri.balance_4tasks  <- df_fmri.balance_4tasks  %>% dplyr::mutate(SN_ECN=SN-ECN, ECN_DMN=ECN-DMN, SN_DMN=SN-DMN) %>% dplyr::filter(!is.na(sex))
df_task <- merge(df_task, df_fmri.balance_4tasks, by.x=c("sub_nr", "Session", "Phase"), by.y=c("sub_nr", "Session", "Run"))

Nback

We start of with the nback task. We first correct the nback trials for the baseline reaction speed using the oddball reaction time from each run as the correcting factor. We then also calculate a measure that combines speed and accuracy after that. Here we use the LIASES method, which is a bit complicated. We do this in a separate chunk.

df_task.nback <- df_task %>% group_by(sub_id, Session, Run) %>% dplyr::mutate(odd_mean_rt=mean(ifelse(Task=="oddball", rt, NA), na.rm=TRUE)) %>% ungroup() %>% subset(Task=="nback") %>% subset(rt>200)
df_task.nback$rt_corr <- df_task.nback$rt/df_task.nback$odd_mean_rt
df_task.nback$Hit[df_task.nback$button_pressed>0 & df_task.nback$button_correct>0] <- 1 
df_task.nback$Miss[df_task.nback$button_pressed==0 & df_task.nback$button_correct>0] <- 1 
df_task.nback$FA[df_task.nback$button_pressed>0 & df_task.nback$button_correct==0] <- 1 
df_task.nback$CorrRej[df_task.nback$button_pressed==0 & df_task.nback$button_correct==0] <- 1 

We now calculate the LIASES score based as detailed in the methods paper by André Vandierendonck (2017). We need to get the RT by participant/condition, the SD of both the RTs and the PEs, and then the total PE per condition. SO we can then compare this measure between runs and conditions. The formula is RTj+ (SD(RT)/SD(PE)) * PEj where j=condition. So I first group by subject, get the RT sd, then group by session and run to get the PE per condition, and the RT per condition, and finally I get the overall PE from all conditions.

df_task.nback <- df_task.nback %>% group_by(sub_id) %>% dplyr::mutate(rt_sd=sd(ifelse(rt!=0, rt, NA), na.rm=TRUE)) %>% group_by(sub_id, Session, Run) %>% dplyr::mutate(pe=((sum(FA, na.rm=TRUE)+ sum(Miss, na.rm=TRUE))/length(rt)), rt_mean=mean(ifelse(rt!=0, rt, NA), na.rm=TRUE)) %>% ungroup %>% group_by(sub_id) %>% dplyr::mutate(pe_sd=sd(pe, na.rm=TRUE)) %>% subset(rt>200) %>% ungroup %>% dplyr::mutate(LIASES=rt_mean+(rt_sd/pe_sd)*pe) %>% ungroup()

Memory

For the memory tasks, we have trials with negative, and others with neutral images. So the goal is to look at RT’s, as well as the percent of correct responses. We subset the memory trials from the main data frame here, and then calculate all the measures we need for analysis.

# Filter out the memory trials
df_task.mem <- df_task %>% subset(Task=="memory")
# Label the factors
df_task.mem$valence <- factor(df_task.mem$valence, levels=c(1,2), labels=c("neutral","negative"))
# Get the number of correct answers and valences for sum calculates
df_task.mem$corr_resp[df_task.mem$button_correct==df_task.mem$button_pressed] <- 1
df_task.mem$incorr_resp[df_task.mem$button_correct!=df_task.mem$button_pressed] <- 1
df_task.mem$mem_neu_n[df_task.mem$valence=="neutral"] <- 1
df_task.mem$mem_neg_n[df_task.mem$valence=="negative"] <- 1
# Get averages and the likes
df_task.mem <- df_task.mem %>% group_by(sub_id, Session, Run) %>% dplyr::mutate(
    mem_neu=sum(ifelse(valence=="neutral", mem_neu_n, NA), na.rm=TRUE),
    mem_neg=sum(ifelse(valence=="negative", mem_neg_n, NA), na.rm=TRUE), 
    mem_tot=mem_neg+mem_neu, 
    mem_neu_corr=sum(ifelse(valence=="neutral", corr_resp, NA), na.rm=TRUE),
    mem_neg_corr=sum(ifelse(valence!="neutral", corr_resp, NA), na.rm=TRUE),
    mem_tot_corr=mem_neg_corr+mem_neu_corr, 
    mem_neg_per=mem_neg_corr/mem_neg,
    mem_neu_per=mem_neu_corr/mem_neu, 
    mem_tot_per=(mem_neu_corr+mem_neg_corr)/mem_tot) %>% ungroup

Oddball

We also need to import the data logs for the oddball. The main outcome measure for this task is the oddball recall done outside the scanner, which we analyze a little differently. For this, we need to calculate dprime. I will also be looking at the RTs here for the individual trials during scanning here. So in addition to the previous data frame from the scan, we also import the data frame from outside the scanner. This was already done in python, so we just need to bring in the new data frame with some minor cleaning. We also calculate dprime for the overall performance, and also the performance by valence

# Oddball df
df_task.odd  <- read.csv("/project/3013068.02/stats/fMRI/Behavioral/fmri_recall.csv", sep="\t", header=TRUE)
#Add variable for separating late, early 
df_task.odd$Scan[df_task.odd$session==1 & df_task.odd$run==1] <- 1
df_task.odd$Scan[df_task.odd$session==1 & df_task.odd$run==2] <-2
df_task.odd$Scan[df_task.odd$session==2 & df_task.odd$run==1] <-3
df_task.odd$Scan[df_task.odd$session==2 & df_task.odd$run==2] <-4
#Make phases into factor
df_task.odd$Session <- factor(df_task.odd$session, levels=c(1,2),
                     labels=c('Control', 'Stress'))
df_task.odd$Scan <- factor(df_task.odd$Scan, levels=c(1,2, 3,4),
                     labels=c('Control-Early','Control-Late',
                              'Stress-Early', 'Stress-Late'))
df_task.odd$Run <- factor(df_task.odd$run, levels=c(1,2),
                     labels=c("Early","Late"))
# Add the MRI signal data
df_task.odd <- merge(df_task.odd, df_fmri.balance_4tasks, by.x=c("sub_nr", "Session", "Run"), by.y=c("sub_nr", "Session", "Run"))
# Now we calculate dprime
##First I drop the NAs
df_task.odd <- na.omit(df_task.odd)
df_task.odd_temp <- df_task.odd
# Then make the frame
df_task.odd_dprime_neg <-as.data.frame(dprime(df_task.odd_temp$odd_hits_neg, df_task.odd_temp$odd_FA_neg, n_miss=df_task.odd_temp$odd_miss_neg, n_cr=df_task.odd_temp$odd_correj_neg))
colnames(df_task.odd_dprime_neg) <- paste(colnames(df_task.odd_dprime_neg), "neg", sep = "_")
df_task.odd <- cbind(df_task.odd, df_task.odd_dprime_neg)
df_task.odd_dprime_pos <- as.data.frame(dprime(df_task.odd_temp$odd_hits_pos, df_task.odd_temp$odd_FA_pos, n_miss=df_task.odd_temp$odd_miss_pos, n_cr=df_task.odd_temp$odd_correj_pos))
colnames(df_task.odd_dprime_pos) <- paste(colnames(df_task.odd_dprime_pos), "pos", sep = "_")
df_task.odd <- cbind(df_task.odd, df_task.odd_dprime_pos)
df_task.odd_dprime <- dprime(df_task.odd_temp$odd_hits, df_task.odd_temp$odd_FA, n_miss=df_task.odd_temp$odd_miss, n_cr=df_task.odd_temp$odd_correj)
df_task.odd <- cbind(df_task.odd, df_task.odd_dprime)

4.2 Stats

# Function to run models in parallel
fx.par_model <- function(signal, dv, iv, re, sub_id, fam, data){
  # Right hand formula; # Left hand formula
  rs_formula=dv; 
  ls_formula= paste(iv, signal,  sep="")
  re_formula=paste("(", re, "|", sub_id,")", sep="")
  # Put together ands run the model
  full_form= paste(rs_formula, paste(ls_formula, re_formula, sep="+"),  sep="~")
  if (fam=="lmer"){
      run_model=lmer(as.formula(full_form),
                     contrasts=list(Session="contr.treatment"),
                     data=data, 
                     control=lmerControl(calc.derivs=FALSE, optCtrl=list(maxfun=100000), optimizer="bobyqa"))
  }else{
  run_model=glmer(as.formula(full_form), 
                  contrasts=list(Session="contr.treatment"),
                  data=data,
                  family=fam,
                  control=glmerControl(calc.derivs=FALSE, optCtrl=list(maxfun=100000), optimizer="bobyqa")) }
  return(run_model)}

Nback

For the nback, we will look at three measures: The reaction time (RT), the number of errors (PE), and the combined score using the LIASES method to combine accuracy and speed. We’ll run the RT model on individual trials where there was a response (so only response trials, without factoring in errors). We then look at the error rate using the proportion of misses and false alarms in the overall number of trials, and finally the LIASES that combines these effects, while also account for within subject deviations. These have all been calculated above, so we can just build the models here, and check the results.

Reaction Time

We first filter out the non-responses, and the responses that are too early.I also z-transform, and rescale the responses to make it easier to apply different types of fits. We then look at the distributions as that will help figure out the optimal models to fit. It appears that there is overall no real effect of reaction time, which I find strange looking at the graphs. Its OK though I guess. We don’t run any post-hoc tests due to no main effects in our model either.

Hist.

# First filter out the 0 RTs
df_task.nback_rts <- subset(df_task.nback, rt>200)
df_task.nback_rts$rt_z <- scale(df_task.nback_rts$rt, center=TRUE, scale=TRUE)
df_task.nback_rts$rt_s <- scales::rescale(df_task.nback_rts$rt_z, to=c(0.01,7))
hist((df_task.nback_rts$rt))

Models

# Contrast to investigate
signals=list("", "ECN")
# Run Models in parallel
models.nback_rt <-  mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="rt_s", 
                                               iv=ifelse(X=="",  "Session*Phase", "1+Session*Phase*"), 
                                               re=ifelse(X=="",  "1+Session*Phase", paste("1+Session*Phase+", X, sep="")), 
                                               sub_id="sub_id", 
                                               fam=Gamma(link="inverse"), 
                                               data=df_task.nback_rts)}, mc.cores=5)

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

#HTML Output
tab.2back_rt <- tab_model(models.nback_rt, dv.labels=c("Signal", "Signal ~ ECN"), show.stat=T, show.se=T)
knitr::asis_output(tab.2back_rt$knitr)

	Signal					Signal ~ ECN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	1.89	0.03	1.83 – 1.95	39.33	<0.001	1.90	0.03	1.83 – 1.96	37.76	<0.001
Session [Stress]	1.00	0.01	0.98 – 1.02	0.26	0.795	1.00	0.01	0.98 – 1.02	-0.07	0.946
Phase1	1.01	0.01	1.00 – 1.02	2.10	0.036	1.01	0.01	0.99 – 1.02	0.99	0.324
Session [Stress] * Phase1	1.00	0.01	0.98 – 1.01	-0.29	0.770	0.99	0.01	0.98 – 1.01	-0.58	0.561
ECN						0.99	0.01	0.97 – 1.00	-1.38	0.166
Session [Stress] * ECN						1.00	0.01	0.98 – 1.02	0.04	0.966
Phase1 * ECN						1.01	0.01	0.99 – 1.02	1.25	0.213
(Session [Stress] Phase1) ECN						1.00	0.01	0.98 – 1.02	0.39	0.694
Random Effects
σ2	0.29					0.29
τ00	0.02 sub_id					0.02 sub_id
τ11	0.00 sub_id.Session1					0.00 sub_id.Session1
	0.00 sub_id.Phase1					0.00 sub_id.Phase1
	0.00 sub_id.Session1:Phase1					0.00 sub_id.ECN
						0.00 sub_id.Session1:Phase1
ρ01	0.07					0.10
	-0.04					-0.05
	0.04					0.09
						0.05
ICC	0.05
N	70 sub_id					70 sub_id
Observations	18456					18456
Marginal R2 / Conditional R2	0.000 / 0.053					0.001 / NA

Post-Hoc

joint_tests(models.nback_rt[[2]], by=c("Session", "Phase"))

Session = Control, Phase = Early:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.046  0.8298

Session = Stress, Phase = Early:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.042  0.8380

Session = Control, Phase = Late:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   3.686  0.0549

Session = Stress, Phase = Late:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   3.444  0.0635

Plot: Main

df_task.nback_summ <- summarySEwithin(data=df_task.nback_rts, measurevar = "rt", idvar="sub_id", withinvars=c("Session", "Phase"), na.rm=FALSE)
df_task.nback_summ$Phase <- as.numeric(df_task.nback_summ$Phase)
df_task.nback_summ$Session <- relevel(df_task.nback_summ$Session, "Stress")
ggplot(df_task.nback_summ, aes(x=Phase, y=rt, color=Session, fill=Session)) +
    geom_line() + 
    geom_errorbar(aes(ymin=rt-se, ymax=rt+se, color=Session),size=1, width=0.20, na.rm=T)+ 
  ggtheme2

Plot: ECN Phase

ggplot(df_task.nback_rts, aes(y=rt_s, x=ECN, colour=Phase)) + geom_smooth(method="lm") + scale_color_manual(values=stress_colors) + ggtheme2

joint_tests(models.nback_rt[[2]],  by=c("Phase"))

Phase = Early:
 model term  df1 df2 F.ratio p.value
 Session       1 Inf   0.112  0.7380
 ECN           1 Inf   0.000  0.9826
 Session:ECN   1 Inf   0.094  0.7596

Phase = Late:
 model term  df1 df2 F.ratio p.value
 Session       1 Inf   0.061  0.8056
 ECN           1 Inf   6.281  0.0122
 Session:ECN   1 Inf   0.052  0.8204

emtrends(models.nback_rt[[2]], pairwise ~ Phase, var="ECN")

NOTE: Results may be misleading due to involvement in interactions

$emtrends
 Phase ECN.trend      SE  df asymp.LCL asymp.UCL
 Early -0.000172 0.00789 Inf   -0.0156   0.01529
 Late  -0.022360 0.00892 Inf   -0.0398  -0.00487

Results are averaged over the levels of: Session 
Confidence level used: 0.95 

$contrasts
 contrast     estimate     SE  df z.ratio p.value
 Early - Late   0.0222 0.0105 Inf   2.106  0.0352

Results are averaged over the levels of: Session 

ggplot(df_task.nback_rts, aes(y=rt_s, x=SN, colour=Session)) + geom_smooth(method="lm") + scale_color_manual(values=stress_colors) + facet_grid(.~Phase) + ggtheme2

Proportion of Errors

For the error rates, we cant really take the trial-by-trial approach like we do with the reaction times, so instead we will take a look at the error rates per run/condition. We can use this filtered data frame later too for LIASES since we have a similar measure per session*run. We first calculate the number of errors instead of the proportion. Thats mainly because the proportion doesn’t fit will in the models, and isn’t an accurate representation given the zero-inflation. So we use counts instead, that way we can use a Poisson model which would fit the data much better.

# Subset
df_task.nback_sum <- df_task.nback_rts
df_task.nback_sum <- df_task.nback_sum %>% dplyr::select("sub_id", "Session", "Phase", "Run","run", "LIASES","pe", "FA", "Miss", "ECN", "SN", "SN_ECN") %>% group_by(sub_id, Session, Run) %>% dplyr::mutate(FA_count=sum(FA, na.rm=TRUE), Miss_count=sum(Miss, na.rm=TRUE), error_count=(FA_count+Miss_count)) %>% ungroup() %>% dplyr::select("sub_id", "Session", "Phase", "Run","run", "LIASES","pe", "error_count", "ECN", "SN", "SN_ECN") %>% unique()
# Rescale the variables
df_task.nback_sum$LIASES_z <- scale(df_task.nback_sum$LIASES, center=TRUE, scale=TRUE)
df_task.nback_sum$pe_z <- scale(df_task.nback_sum$pe, center=TRUE, scale=TRUE)

Models

# Contrasts to investigate
signals=list("", "ECN")
# Ru moddels in parallel
models.nback_pe <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="error_count", 
                                               iv=ifelse(X=="",  "Session*Phase", "1+Session*Phase*"), 
                                               re=ifelse(X=="",  "1+Session*Phase", paste("1+Session*Phase+", X, sep="")), 
                                               sub_id="sub_id", 
                                               fam=poisson(), 
                                               data=df_task.nback_sum)}, mc.cores=5 )

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

#HTML Output
tab.2back_pe <- tab_model(models.nback_pe, dv.labels=signals, show.se=T, show.stat=T)
asis_output(tab.2back_pe$knitr)

						ECN
Predictors	Incidence Rate Ratios	std. Error	CI	Statistic	p	Incidence Rate Ratios	std. Error	CI	Statistic	p
(Intercept)	0.99	0.10	0.82 – 1.20	-0.09	0.931	0.95	0.10	0.77 – 1.18	-0.46	0.646
Session [Stress]	0.89	0.07	0.76 – 1.04	-1.51	0.131	0.95	0.09	0.79 – 1.15	-0.53	0.595
Phase1	1.01	0.05	0.92 – 1.10	0.16	0.875	0.99	0.06	0.87 – 1.11	-0.24	0.812
Session [Stress] * Phase1	0.96	0.07	0.84 – 1.11	-0.53	0.594	1.01	0.09	0.85 – 1.19	0.12	0.908
ECN						1.06	0.09	0.91 – 1.25	0.76	0.444
Session [Stress] * ECN						0.91	0.09	0.74 – 1.11	-0.95	0.345
Phase1 * ECN						1.03	0.07	0.90 – 1.18	0.46	0.647
(Session [Stress] Phase1) ECN						0.96	0.09	0.80 – 1.16	-0.43	0.668
Random Effects
σ2	0.72					0.72
τ00	0.58 sub_id					0.57 sub_id
τ11	0.03 sub_id.Session1					0.03 sub_id.Session1
	0.00 sub_id.Phase1					0.01 sub_id.Phase1
	0.01 sub_id.Session1:Phase1					0.03 sub_id.ECN
						0.01 sub_id.Session1:Phase1
ρ01	-0.46					-0.47
	-0.08					-0.23
	0.20					-0.17
						0.41
ICC						0.44
N	70 sub_id					70 sub_id
Observations	838					838
Marginal R2 / Conditional R2	0.005 / NA					0.004 / 0.447

Post-Hoc: Main

emmeans(models.nback_pe[[1]],  list(pairwise ~Session*Phase), by=c("Phase"), adjust = "tukey")

$`emmeans of Session | Phase`
Phase = Early:
 Session   emmean    SE  df asymp.LCL asymp.UCL
 Control -0.00137 0.110 Inf    -0.217    0.2140
 Stress  -0.15670 0.124 Inf    -0.399    0.0858

Phase = Late:
 Session   emmean    SE  df asymp.LCL asymp.UCL
 Control -0.01545 0.103 Inf    -0.217    0.1859
 Stress  -0.09563 0.127 Inf    -0.344    0.1527

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of Session | Phase`
Phase = Early:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.1553 0.1103 Inf   1.408  0.1590

Phase = Late:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.0802 0.0996 Inf   0.805  0.4208

Results are given on the log (not the response) scale. 

Post-Hoc: Brain

joint_tests(models.nback_pe[[2]], by=c("Session", "Phase")) 

Session = Control, Phase = Early:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.795  0.3725

Session = Stress, Phase = Early:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.199  0.6554

Session = Control, Phase = Late:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.082  0.7745

Session = Stress, Phase = Late:
 model term df1 df2 F.ratio p.value
 ECN          1 Inf   0.052  0.8199

Plot

df.nback_error_sum <- summarySEwithin(data=df_task.nback_sum, measurevar="error_count", idvar="sub_id", withinvars=  c("Phase", "Session"), na.rm=FALSE)
df.nback_error_sum$Phase <- as.numeric(df.nback_error_sum$Phase)
df.nback_error_sum$Session <- relevel(df.nback_error_sum$Session, "Stress")
plot.nback_errors <- ggplot(df.nback_error_sum, aes(x=Phase, y=error_count, color=Session, fill=Session)) +
    geom_bar(stat="summary", position="dodge2", alpha=0.5) + 
    geom_errorbar(aes(ymin=error_count-se, ymax=error_count+se, color=Session),size=1, width=0.45, na.rm=T, position=position_dodge(.9) ) + ggtheme2
plot.nback_errors

LIASES

The final measure we use to calculate the n-back performance is LIASES, which combines accuracy and speed. This method is mentioned earlier, and a data frame was subsetted in the chunk looking at the error rates, so we just need to build the model and run it here. I first tested different fits, and it appears that the inverse Gaussian is the best to fit this data.

Hist.

hist((df_task.nback_sum$LIASES_z))

Models

# Parallel
signals=list("",  "ECN" )
models.nback_liases <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="LIASES_z",  
                                               iv=ifelse(X=="",  "Session*Phase", "1+Session*Phase*"), 
                                               re=ifelse(X=="",  "1+Session*Phase", paste("1+Session*Phase+", X, sep="")), 
                                               sub_id="sub_id", 
                                               fam=gaussian(link="inverse"), 
                                               data=df_task.nback_sum)}, mc.cores=5)
#HTML Output
tab.2back_liases <- tab_model(models.nback_liases, dv.labels=c("Signal","Signal ~ ECN" ),show.stat=T, show.se=T)
asis_output(tab.2back_liases$knitr)

	Signal					Signal ~ ECN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.17	0.21	-0.25 – 0.59	0.82	0.415	0.28	0.27	-0.24 – 0.81	1.05	0.292
Session [Stress]	-0.29	0.21	-0.71 – 0.13	-1.36	0.174	-0.22	0.26	-0.73 – 0.29	-0.84	0.399
Phase1	0.02	0.13	-0.24 – 0.28	0.12	0.904	-0.09	0.17	-0.42 – 0.24	-0.56	0.579
Session [Stress] * Phase1	-0.20	0.21	-0.61 – 0.21	-0.96	0.338	-0.02	0.24	-0.48 – 0.44	-0.08	0.935
ECN						-0.14	0.19	-0.52 – 0.23	-0.75	0.455
Session [Stress] * ECN						-0.31	0.26	-0.82 – 0.21	-1.17	0.240
Phase1 * ECN						0.15	0.17	-0.18 – 0.48	0.90	0.368
(Session [Stress] Phase1) ECN						-0.16	0.25	-0.64 – 0.33	-0.63	0.528
Random Effects
σ2	0.51					0.51
τ00	2.13 sub_id					2.57 sub_id
τ11	0.57 sub_id.Session1					0.53 sub_id.Session1
	0.44 sub_id.Phase1					0.46 sub_id.Phase1
	0.53 sub_id.Session1:Phase1					0.12 sub_id.ECN
						0.38 sub_id.Session1:Phase1
ρ01	0.02					0.09
	-0.12					-0.12
	0.12					-0.74
						0.17
ICC	0.83					0.84
N	70 sub_id					70 sub_id
Observations	838					838
Marginal R2 / Conditional R2	0.012 / 0.836					0.039 / 0.845

Post-Hoc: Behavioural

emmeans(models.nback_liases[[1]], pairwise ~ Session | Phase, adjust="none")

Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates

$emmeans
Phase = Early:
 Session  emmean    SE  df asymp.LCL asymp.UCL
 Control  0.1907 0.264 Inf    -0.327     0.709
 Stress  -0.2990 0.232 Inf    -0.754     0.156

Phase = Late:
 Session  emmean    SE  df asymp.LCL asymp.UCL
 Control  0.1586 0.239 Inf    -0.311     0.628
 Stress   0.0675 0.286 Inf    -0.492     0.627

Results are given on the inverse (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Phase = Early:
 contrast         estimate    SE  df z.ratio p.value
 Control - Stress    0.490 0.299 Inf   1.639  0.1012

Phase = Late:
 contrast         estimate    SE  df z.ratio p.value
 Control - Stress    0.091 0.297 Inf   0.306  0.7595

Note: contrasts are still on the inverse scale 

Plot: Behavioral

df_task.nback_summ <- summarySEwithin(df_task.nback, measurevar="LIASES", idvar="sub_id", withinvar=c("Session", "Phase"), na.rm=T)

Automatically converting the following non-factors to factors: Phase

df_task.nback_summ$run <- as.numeric(df_task.nback_summ$Phase)
df_task.nback_summ$Session <- relevel(df_task.nback_summ$Session, "Stress")
ggplot(df_task.nback_summ, aes(x=run, y=LIASES, color=Session, fill=Session)) +
    geom_line() + 
    geom_errorbar(aes(ymin=LIASES-se, ymax=LIASES+se, color=Session),size=0.5, width=0.15, na.rm=T) + 
  ggtheme2 

Post-Hoc: ECN by Session

emtrends(models.nback_liases[[2]], identity ~ Session, var="ECN")

NOTE: Results may be misleading due to involvement in interactions

$emtrends
 Session ECN.trend    SE  df asymp.LCL asymp.UCL
 Control    -0.144 0.192 Inf     -0.52    0.2330
 Stress     -0.451 0.214 Inf     -0.87   -0.0318

Results are averaged over the levels of: Phase 
Confidence level used: 0.95 

$contrasts
 contrast estimate    SE  df z.ratio p.value
 Control    -0.144 0.192 Inf  -0.748  0.4547
 Stress     -0.451 0.214 Inf  -2.109  0.0350

Results are averaged over the levels of: Phase 

#emtrends(models.nback_liases[[3]], pairwise ~ Session | Phase, var="ECN")

ggplot(df_task.nback_sum, aes(x=ECN, y=LIASES_z, colour=Session)) +
  geom_smooth(method="lm", se=T) + 
  xlab("ECN Activity") + ylab("LIASES Accuracy Score") +
  ggtheme2  + scale_color_manual(values=stress_colors)

Memory

For the memory task, we look at the reaction time to the stimuli, split by valence, and the recognition of the images we presented them with outside the scanner. For the RT data, we will take a look at individual trials, and then for the recognition of the images we take a look at the average reposes per run or phase, depending on the model fits best. We first also remove any rts that are less than 200ms as those are just too quick, and we also rescale the RT variables.

df_task.mem <- df_task.mem %>% subset(rt>=200)
df_task.mem$rt_z <- scale(df_task.mem$rt)
df_task.mem$rt_s <- scales::rescale(df_task.mem$rt_z, to=c(1,6))

Recognition

We next look at the ability to recall the object-location associations. We first will look at overall recall, then take a look at the recall by valence associations to see if the negative or positively valent images are remembered better. First I subset the data, and make the data frame a bit longer so that we have percent correct for each of the neutral and the negative items in one column.

# Subset the neutral item correct response rates
df_task.mem_neu <- df_task.mem %>% dplyr::select("sub_id", "Session", "Phase", "Run","run", "mem_neu_per", "DMN") %>% unique() %>% dplyr::rename(per_cor=mem_neu_per) %>% mutate(valence="neu")
# Do the same for the negative items
df_task.mem_neg <- df_task.mem %>% dplyr::select("sub_id", "Session", "Phase", "Run","run", "mem_neg_per", "DMN") %>% unique() %>% dplyr::rename(per_cor=mem_neg_per)  %>% mutate(valence="neg")
# Bind them into one dataframe
df_task.mem_cor <- rbind.data.frame(df_task.mem_neg, df_task.mem_neu)

Now we can also check out the distribution of the correct response rates for each of the valences to help pick the most suitable model family. We see that both valences look similar, but skewed left.

Hist.

hist(df_task.mem_neu$per_cor)

hist(df_task.mem_neg$per_cor)

Chance Level

model.memory_chance <- t.test(df_task.mem$mem_tot_per, mu=0.25)
tidy(model.memory_chance)

Models

# Main
model.mem_correct <- lmer(per_cor ~ Session*Phase*valence + (1+Session*Phase|sub_id),
                          contrasts=list(Session="contr.treatment"),
                           data=df_task.mem_cor,
                           control=lmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)) )
#check_model(model.mem_correct)
# Brain
model.mem_correct_brain <- lmer(per_cor ~ Session*Phase*DMN + valence + (1+Session*Phase*DMN|sub_id),
                                contrasts=list(Session="contr.treatment"),
                           data=df_task.mem_cor,
                           control=lmerControl(calc.derivs = FALSE, optimizer="bobyqa", optCtrl=list(maxfun=100000)) )
# HTM Table
asis_output(tab_model(model.mem_correct, model.mem_correct_brain)$knitr)

	per cor			per cor
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.61	0.57 – 0.64	<0.001	0.62	0.58 – 0.66	<0.001
Session [Stress]	0.01	-0.01 – 0.03	0.474	0.00	-0.02 – 0.03	0.841
Phase1	-0.01	-0.02 – 0.01	0.581	-0.01	-0.03 – 0.01	0.510
valence1	0.03	0.02 – 0.04	<0.001	0.03	0.02 – 0.03	<0.001
Session [Stress] * Phase1	0.01	-0.01 – 0.03	0.306	0.01	-0.02 – 0.04	0.421
Session [Stress] * valence1	-0.00	-0.02 – 0.01	0.856
Phase1 * valence1	-0.00	-0.01 – 0.01	0.440
Session [Stress] * Phase1 * valence1	0.00	-0.01 – 0.02	0.853
DMN				0.02	-0.00 – 0.05	0.099
Session [Stress] * DMN				-0.01	-0.04 – 0.02	0.575
Phase1 * DMN				-0.01	-0.04 – 0.02	0.623
(Session [Stress] Phase1) DMN				0.00	-0.03 – 0.04	0.827
Random Effects
σ2	0.02			0.02
τ00	0.02 sub_id			0.02 sub_id
τ11	0.00 sub_id.Session1			0.00 sub_id.Session1
	0.00 sub_id.Phase1			0.00 sub_id.Phase1
	0.00 sub_id.Session1:Phase1			0.00 sub_id.DMN
				0.00 sub_id.Session1:Phase1
				0.00 sub_id.Session1:DMN
				0.00 sub_id.Phase1:DMN
				0.00 sub_id.Session1:Phase1:DMN
ρ01	-0.04			0.04
	-0.31			-0.45
	0.04			0.61
				-0.09
				0.63
				-0.37
				-0.29
ICC	0.48
N	70 sub_id			70 sub_id
Observations	1676			1676
Marginal R2 / Conditional R2	0.015 / 0.489			0.033 / NA

Post-Hoc: Main

Only the valence effect sticks out, so we run a specific post-hoc test to investigate it further.

emmeans(model.mem_correct,  list(pairwise ~Session*valence), by=c("Phase"), adjust = "tukey")

$`emmeans of Session, valence | Phase`
Phase = Early:
 Session valence emmean     SE   df lower.CL upper.CL
 Control neg      0.625 0.0215 90.5    0.582    0.667
 Stress  neg      0.645 0.0206 92.8    0.604    0.686
 Control neu      0.581 0.0215 90.5    0.538    0.624
 Stress  neu      0.601 0.0206 92.8    0.561    0.642

Phase = Late:
 Session valence emmean     SE   df lower.CL upper.CL
 Control neg      0.644 0.0230 87.3    0.598    0.689
 Stress  neg      0.637 0.0228 87.5    0.592    0.683
 Control neu      0.583 0.0230 87.3    0.538    0.629
 Stress  neu      0.582 0.0228 87.5    0.537    0.628

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of Session, valence | Phase`
Phase = Early:
 2                          estimate     SE   df t.ratio p.value
 Control neg - Stress neg  -0.020473 0.0211  126  -0.969  0.7671
 Control neg - Control neu  0.043547 0.0153 1392   2.839  0.0237
 Control neg - Stress neu   0.023124 0.0211  126   1.095  0.6934
 Stress neg - Control neu   0.064020 0.0211  126   3.031  0.0155
 Stress neg - Stress neu    0.043597 0.0153 1392   2.849  0.0230
 Control neu - Stress neu  -0.020423 0.0211  126  -0.967  0.7685

Phase = Late:
 2                          estimate     SE   df t.ratio p.value
 Control neg - Stress neg   0.006374 0.0185  158   0.345  0.9858
 Control neg - Control neu  0.060314 0.0153 1392   3.932  0.0005
 Control neg - Stress neu   0.061059 0.0185  158   3.305  0.0064
 Stress neg - Control neu   0.053940 0.0185  158   2.920  0.0207
 Stress neg - Stress neu    0.054685 0.0153 1392   3.574  0.0021
 Control neu - Stress neu   0.000746 0.0185  158   0.040  1.0000

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates 

Plot: Main

df_task.mem_cor_sum <- summarySEwithin(data=df_task.mem_cor, idvar="sub_id", measurevar="per_cor", withinvars=c("Session","Phase","valence"))

Automatically converting the following non-factors to factors: Phase, valence

plot.mem_corr <- ggplot(df_task.mem_cor_sum, aes(y=per_cor, x=Phase, colour=valence, fill=valence, na.rm = TRUE)) +
  geom_bar(stat="summary", alpha=.5, position=position_dodge2())+
  geom_errorbar(aes(ymin=per_cor-se, ymax=per_cor+se), position=position_dodge(.9),size=1, width=0.45, na.rm=T) +
  scale_y_continuous()+ scale_x_discrete()+
  ggtitle("Picture-Location Recall")+
  ylab("Percent Remembered\n") + 
  ggtheme2 
plot.mem_corr + facet_grid(.~Session)

Post-Hoc: DMN by Phase

emtrends(model.mem_correct_brain, pairwise ~ Phase, var="DMN") 

NOTE: Results may be misleading due to involvement in interactions

$emtrends
 Phase DMN.trend     SE   df lower.CL upper.CL
 Early    0.0112 0.0177 28.0  -0.0250   0.0474
 Late     0.0213 0.0178 28.9  -0.0152   0.0578

Results are averaged over the levels of: Session, valence 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast     estimate     SE   df t.ratio p.value
 Early - Late  -0.0101 0.0252 33.9  -0.399  0.6928

Results are averaged over the levels of: Session, valence 
Degrees-of-freedom method: kenward-roger 

emmip(model.mem_correct_brain, Phase ~ DMN, cov.reduce = range) + scale_color_manual(values=stress_colors) + ggtheme2

NOTE: Results may be misleading due to involvement in interactions

Reaction Time

We first take a look at the reaction time. This analysis is done on a trial by trial basis. We will look at the RT between the stress and control sessions, and also look at the effects of the valence of the images on the reaction time. First we take a look at the distribution of the RT’s to better inform us of which model family would fit best. It appears here that we have a significant effect of valence on reaction times. So participants are faster in responding.

Hist.

hist(df_task.mem$rt)

Models

# Main
model.memory_rt <- glmer(rt_s ~ Session*Phase + Session*valence +Session*Phase*valence + (1+Session*Phase|sub_id),
                         contrasts=list(Session="contr.treatment"),
                        data=df_task.mem,
                        family=Gamma("inverse"),
                       control=glmerControl(calc.derivs=FALSE))
# Brain
model.memory_rt_brain <- glmer(rt_s ~ Session*Phase*DMN + (1+Session*Phase*DMN|sub_id),
                               contrasts=list(Session="contr.treatment"),
                        data=df_task.mem,
                       family=Gamma("identity"),
                       control=glmerControl(calc.derivs=FALSE, optimizer="bobyqa", optCtrl=list(maxfun=100000)))
# HTML Output
asis_output(tab_model(model.memory_rt, model.memory_rt_brain)$knitr)

	rt s			rt s
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	1.62	1.60 – 1.63	<0.001	2.19	2.16 – 2.23	<0.001
Session [Stress]	1.00	1.00 – 1.01	0.474	-0.01	-0.04 – 0.03	0.753
Phase1	1.00	0.99 – 1.00	0.038	0.03	0.01 – 0.05	0.001
valence1	1.00	1.00 – 1.00	0.062
Session [Stress] * Phase1	1.00	1.00 – 1.01	0.438	-0.02	-0.05 – 0.01	0.220
Session [Stress] * valence1	1.00	1.00 – 1.01	0.368
Phase1 * valence1	1.00	1.00 – 1.00	0.195
Session [Stress] * Phase1 * valence1	1.00	1.00 – 1.00	0.935
DMN				-0.02	-0.05 – 0.01	0.191
Session [Stress] * DMN				0.03	-0.01 – 0.08	0.146
Phase1 * DMN				0.05	0.02 – 0.08	0.003
(Session [Stress] Phase1) DMN				-0.01	-0.05 – 0.04	0.703
Random Effects
σ2	0.07			0.07
τ00	0.00 sub_id			0.01 sub_id
τ11	0.00 sub_id.Session1			0.00 sub_id.Session1
	0.00 sub_id.Phase1			0.00 sub_id.Phase1
	0.00 sub_id.Session1:Phase1			0.00 sub_id.DMN
				0.00 sub_id.Session1:Phase1
				0.00 sub_id.Session1:DMN
				0.00 sub_id.Phase1:DMN
				0.00 sub_id.Session1:Phase1:DMN
ρ01	-0.01			-0.01
	-0.04			-0.06
	0.11			-0.96
				0.15
				-0.12
				0.24
				0.23
ICC	0.01
N	70 sub_id			70 sub_id
Observations	17529			17529
Marginal R2 / Conditional R2	0.000 / 0.011			0.020 / NA

Post-Hoc: Main

emmeans(model.memory_rt,  list(pairwise ~ Session*Phase), by=c("Phase"), adjust = "tukey")

NOTE: Results may be misleading due to involvement in interactions
Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates

$`emmeans of Session | Phase`
Phase = Early:
 Session emmean      SE  df asymp.LCL asymp.UCL
 Control  0.475 0.00417 Inf     0.467     0.483
 Stress   0.480 0.00409 Inf     0.472     0.488

Phase = Late:
 Session emmean      SE  df asymp.LCL asymp.UCL
 Control  0.484 0.00425 Inf     0.476     0.493
 Stress   0.485 0.00436 Inf     0.476     0.493

Results are averaged over the levels of: valence 
Results are given on the inverse (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of Session | Phase`
Phase = Early:
 2                 estimate      SE  df z.ratio p.value
 Control - Stress -0.004940 0.00455 Inf  -1.085  0.2780

Phase = Late:
 2                 estimate      SE  df z.ratio p.value
 Control - Stress -0.000372 0.00491 Inf  -0.076  0.9396

Results are averaged over the levels of: valence 
Note: contrasts are still on the inverse scale 

Plot: Main

df_task.mem_sum <- summarySEwithin(df_task.mem, measurevar="rt",idvar="sub_id", withinvars= c("Session","Phase","valence"))

Automatically converting the following non-factors to factors: Phase

df_task.mem_sum$run = as.numeric((df_task.mem_sum$Phase))
df_task.mem_sum$valence <- relevel(df_task.mem_sum$valence, "negative")
plot.mem_rt <- ggplot(df_task.mem_sum, aes(x=run, y=rt, colour=valence)) + 
                          geom_line() +
                          geom_errorbar(aes(ymin=rt+se, ymax=rt-se),width=0.2, na.rm=TRUE) + 
                          scale_x_continuous(breaks=c(1,2), labels = c("Early", "Late")) + 
                          ggtheme2  
plot.mem_rt + facet_grid(.~Session)

Post-Hoc: Brain

joint_tests(model.memory_rt_brain, by=c("Session","Phase"))

Session = Control, Phase = Early:
 model term df1 df2 F.ratio p.value
 DMN          1 Inf   1.768  0.1836

Session = Stress, Phase = Early:
 model term df1 df2 F.ratio p.value
 DMN          1 Inf   4.950  0.0261

Session = Control, Phase = Late:
 model term df1 df2 F.ratio p.value
 DMN          1 Inf   7.306  0.0069

Session = Stress, Phase = Late:
 model term df1 df2 F.ratio p.value
 DMN          1 Inf   1.464  0.2264

ggplot(df_task.mem, aes(x=DMN, y=rt_s, colour=Session)) + geom_smooth(method="lm") + ggtheme2 +  facet_grid(.~Phase) + scale_color_manual(values=stress_colors)

Oddball

The oddball data has already been processed, theres not much more we need to do there. Therefore, we can begin with the stats.

d-Prime

Hist.

hist(df_task.odd$dprime)

Models

# Signals
signals=c("", "SN" )
# Run in parallel function
models.odd_dprime <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="dprime",  
                                               iv=ifelse(X=="",  "Session*Run", "1+Session*Run*"), 
                                               re=ifelse(X=="",  "1+Session+Run", paste("1+Session+Run")), 
                                               sub_id="sub_nr", 
                                               fam="lmer", 
                                               data=df_task.odd)}, mc.cores=5)
# Print HTML Output
asis_output(tab_model(models.odd_dprime, dv.labels = c("Signal", "Signal ~ SN"), show.stat =T, show.se=T)$knitr)

	Signal					Signal ~ SN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.40	0.05	0.31 – 0.49	8.58	<0.001	0.34	0.07	0.20 – 0.47	4.92	<0.001
Session [Stress]	-0.43	0.06	-0.54 – -0.32	-7.39	<0.001	-0.38	0.09	-0.56 – -0.20	-4.10	<0.001
Run1	0.01	0.03	-0.05 – 0.07	0.38	0.703	0.04	0.06	-0.07 – 0.15	0.69	0.492
Session [Stress] * Run1	-0.04	0.04	-0.12 – 0.05	-0.90	0.369	-0.08	0.08	-0.24 – 0.07	-1.04	0.299
SN						0.05	0.05	-0.04 – 0.15	1.15	0.250
Session [Stress] * SN						-0.05	0.06	-0.17 – 0.08	-0.75	0.452
Run1 * SN						-0.02	0.04	-0.10 – 0.06	-0.43	0.667
(Session [Stress] * Run1) * SN						0.03	0.06	-0.08 – 0.15	0.56	0.575
Random Effects
σ2	0.10					0.10
τ00	0.02 sub_nr					0.02 sub_nr
τ11	0.02 sub_nr.Session1					0.02 sub_nr.Session1
	0.00 sub_nr.Run1					0.00 sub_nr.Run1
ρ01	0.57					0.58
	-0.82					-0.88
N	63 sub_nr					63 sub_nr
Observations	222					222
Marginal R2 / Conditional R2	0.319 / NA					0.317 / NA

Plot: Main

plot.odd_sess <- ggplot(df_task.odd, aes(x=Scan, y=dprime, colour=Scan, fill=Scan)) + 
    geom_violin(alpha=0.5) + 
    geom_boxplot(width=0.1, alpha=0.2) +
    geom_hline(yintercept=0, colour="grey") +
    ggtheme2
plot.odd_sess

Post-Hoc: Brain

joint_tests(models.odd_dprime[[2]], by=c("Session", "Run"))

Session = Control, Run = Early:
 model term df1    df2 F.ratio p.value
 SN           1  87.57   0.328  0.5686

Session = Stress, Run = Early:
 model term df1    df2 F.ratio p.value
 SN           1  85.42   0.118  0.7318

Session = Control, Run = Late:
 model term df1    df2 F.ratio p.value
 SN           1 101.58   1.121  0.2922

Session = Stress, Run = Late:
 model term df1    df2 F.ratio p.value
 SN           1  86.14   0.022  0.8818

Hits

Hist.

We first check the hits. It appears that there is a significant difference in the hit-rate between the two sessions, and between the runs.

hist(df_task.odd$odd_hits)

Models

# Contrasts to investigate
signals=c("", "SN")
# Run in parallel function
models.odd_hits <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="odd_hits",  
                                               iv=ifelse(X=="",  "Session*Run", "1+Session*Run*"), 
                                               re=ifelse(X=="",  "1+Session+Run", paste("1+Session+Run")), 
                                               sub_id="sub_nr", 
                                               fam="lmer", 
                                               data=df_task.odd)}, mc.cores=5)
# Print HTML Output
asis_output(tab_model(models.odd_hits, dv.labels = c("Signal", "Signal ~ SN"), show.stat =T, show.se=T)$knitr)

	Signal					Signal ~ SN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	20.34	0.69	18.99 – 21.69	29.46	<0.001	20.21	0.95	18.35 – 22.07	21.34	<0.001
Session [Stress]	-3.40	0.68	-4.74 – -2.06	-4.97	<0.001	-3.88	1.26	-6.36 – -1.41	-3.08	0.002
Run1	-1.85	0.43	-2.68 – -1.01	-4.34	<0.001	-1.77	0.73	-3.20 – -0.33	-2.42	0.016
Session [Stress] * Run1	0.22	0.57	-0.90 – 1.34	0.38	0.701	1.06	1.11	-1.12 – 3.24	0.95	0.340
SN						0.12	0.59	-1.04 – 1.28	0.20	0.842
Session [Stress] * SN						0.50	0.92	-1.30 – 2.30	0.54	0.587
Run1 * SN						-0.07	0.55	-1.14 – 1.00	-0.13	0.900
(Session [Stress] * Run1) * SN						-0.79	0.83	-2.42 – 0.83	-0.96	0.339
Random Effects
σ2	17.02					17.40
τ00	24.18 sub_nr					24.36 sub_nr
τ11	2.15 sub_nr.Session1					2.17 sub_nr.Session1
	1.10 sub_nr.Run1					0.87 sub_nr.Run1
ρ01	-0.46					-0.45
	-0.11					-0.11
ICC	0.60					0.59
N	63 sub_nr					63 sub_nr
Observations	222					222
Marginal R2 / Conditional R2	0.126 / 0.649					0.132 / 0.646

Plot: Main

plot.odd_hits <- ggplot(df_task.odd, aes(x=Scan, y=odd_hits, colour=Scan, fill=Scan)) + 
    geom_violin(alpha=0.5) + 
    geom_boxplot(width=0.1, alpha=0.2) +
    ggtheme2
plot.odd_hits 

Post-Hoc: Brain

joint_tests(models.odd_hits[[2]], by=c("Session", "Run")) 

Session = Control, Run = Early:
 model term df1   df2 F.ratio p.value
 SN           1 88.29   0.004  0.9526

Session = Stress, Run = Early:
 model term df1   df2 F.ratio p.value
 SN           1 97.52   0.063  0.8024

Session = Control, Run = Late:
 model term df1   df2 F.ratio p.value
 SN           1 97.12   0.048  0.8266

Session = Stress, Run = Late:
 model term df1   df2 F.ratio p.value
 SN           1 98.67   2.142  0.1465

False Alarms

We also see that they perform more errors. So we have decreased hits, and more errors in the stress session

Models

# MLoop over
signals=c("", "SN")
# Run in parallel function
models.odd_fas <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="odd_FA",  
                                               iv=ifelse(X=="",  "Session*Run", "1+Session*Run*"), 
                                               re=ifelse(X=="",  "1+Session+Run", paste("1+Session+Run")), 
                                               sub_id="sub_nr", 
                                               fam="lmer", 
                                               data=df_task.odd)}, mc.cores=5)
# HTML Output
tab.odd_fas <- tab_model(models.odd_fas, show.se=T, show.stat=T, dv.labels=c("Signal", "Signal ~ SN"))
asis_output(tab.odd_fas$knitr)

	Signal					Signal ~ SN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	7.40	0.36	6.70 – 8.11	20.60	<0.001	8.08	0.50	7.09 – 9.07	16.03	<0.001
Session [Stress]	1.67	0.43	0.82 – 2.52	3.86	<0.001	0.62	0.72	-0.79 – 2.03	0.86	0.391
Run1	-1.07	0.23	-1.52 – -0.62	-4.65	<0.001	-1.45	0.41	-2.25 – -0.65	-3.54	<0.001
Session [Stress] * Run1	0.31	0.31	-0.30 – 0.93	1.01	0.315	1.17	0.62	-0.05 – 2.38	1.88	0.060
SN						-0.58	0.34	-1.24 – 0.08	-1.73	0.083
Session [Stress] * SN						0.95	0.51	-0.05 – 1.95	1.85	0.064
Run1 * SN						0.28	0.30	-0.31 – 0.88	0.93	0.352
(Session [Stress] * Run1) * SN						-0.74	0.46	-1.64 – 0.16	-1.61	0.107
Random Effects
σ2	5.12					5.54
τ00	4.34 sub_nr					4.12 sub_nr
τ11	1.38 sub_nr.Session1					1.07 sub_nr.Session1
	0.23 sub_nr.Run1					0.13 sub_nr.Run1
ρ01	-0.15					-0.24
	-0.21					-0.20
ICC	0.47					0.44
N	63 sub_nr					63 sub_nr
Observations	222					222
Marginal R2 / Conditional R2	0.135 / 0.544					0.151 / 0.521

Plot: Main

plot.odd_fas <- ggplot(df_task.odd, aes(x=Scan, y=odd_FA, colour=Scan, fill=Scan)) + 
    geom_violin(alpha=0.5) + 
    geom_boxplot(width=0.1, alpha=0.2) +
    ggtheme2
plot.odd_fas

Post-Hoc: SN by Session

emtrends(models.odd_fas[[2]], pairwise ~ Session, var="SN")

NOTE: Results may be misleading due to involvement in interactions

$emtrends
 Session SN.trend    SE  df lower.CL upper.CL
 Control   -0.582 0.350 142   -1.275     0.11
 Stress     0.366 0.406 149   -0.436     1.17

Results are averaged over the levels of: Run 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast         estimate    SE  df t.ratio p.value
 Control - Stress   -0.948 0.532 193  -1.783  0.0762

Results are averaged over the levels of: Run 
Degrees-of-freedom method: kenward-roger 

joint_tests(models.odd_fas[[2]], by=c("Session"))

Session = Control:
 model term df1    df2 F.ratio p.value
 Run          1 111.50  21.681  <.0001
 SN           1 142.00   2.761  0.0988
 Run:SN       1 144.07   0.800  0.3727

Session = Stress:
 model term df1    df2 F.ratio p.value
 Run          1 111.64  11.243  0.0011
 SN           1 149.25   0.813  0.3688
 Run:SN       1 139.00   1.636  0.2030

plot.sn_fas <- ggplot(df_task.odd, aes(x=SN, y=odd_FA, colour=Session, fill=Session)) + 
  geom_smooth(method="lm", formula="y~x", alpha=0.15) + 
  #geom_point(alpha=0.5)+
  scale_color_manual(values=stress_colors) +  scale_fill_manual(values=stress_colors) + 
  ylab("False Alarms") +
  ggtheme2
plot.sn_fas + facet_grid(~Run)

dPrime by Valence

Now we take a look at the effects of stress and valence on recall. So we do something similar here to what we did before, where we split the positive and negative values, and then merge them into a long data frame(same as we did for the memory task recognition trials).

# Subset the neutral item correct response rates
df_task.odd_pos <- df_task.odd %>% dplyr::select("sub_nr", "Session", "Run", "Scan","dprime_pos", "SN", "ECN", "SN_ECN") %>% unique() %>% dplyr::rename(dprime=dprime_pos) %>% mutate(Valence="Positive")
# Do the same for the negative items
df_task.odd_neg <- df_task.odd %>% dplyr::select("sub_nr", "Session", "Run", "Scan","dprime_neg", "SN", "ECN", "SN_ECN") %>% unique() %>% dplyr::rename(dprime=dprime_neg) %>% mutate(Valence="Negative")
# Bind them into one dataframe
df_task.odd_val <- rbind.data.frame(df_task.odd_pos, df_task.odd_neg)

Models

# Loop over these
signals=c("", "SN")
# Run in parallel function
models.odd_val <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X, 
                                               dv="dprime",  
                                               iv=ifelse(X=="",  "Session*Run*Valence", "1+Session*Run*Valence*"), 
                                               re=ifelse(X=="",  "1+Session+Run", paste("1+Session+Run")), 
                                               sub_id="sub_nr", 
                                               fam="lmer", 
                                               data=df_task.odd_val)}, mc.cores=5)
# Print HTML Output
asis_output(tab_model(models.odd_val, 
                      dv.labels = c("Signal", "Signal ~ SN"),
                      show.stat =T, show.se=T)$knitr)

	Signal					Signal ~ SN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	0.39	0.05	0.30 – 0.48	8.49	<0.001	0.33	0.07	0.20 – 0.46	4.86	<0.001
Session [Stress]	-0.42	0.06	-0.54 – -0.31	-7.37	<0.001	-0.38	0.09	-0.56 – -0.20	-4.16	<0.001
Run1	0.01	0.03	-0.05 – 0.07	0.34	0.736	0.04	0.06	-0.07 – 0.14	0.67	0.505
Valence1	-0.01	0.03	-0.07 – 0.04	-0.48	0.633	0.01	0.05	-0.09 – 0.10	0.12	0.904
Session [Stress] * Run1	-0.03	0.04	-0.12 – 0.05	-0.76	0.449	-0.08	0.08	-0.24 – 0.08	-0.98	0.327
Session [Stress] * Valence1	0.04	0.04	-0.04 – 0.12	0.91	0.363	0.05	0.07	-0.09 – 0.20	0.73	0.462
Run1 * Valence1	0.03	0.03	-0.02 – 0.09	1.15	0.249	0.01	0.05	-0.09 – 0.11	0.24	0.808
Session [Stress] * Run1 * Valence1	-0.02	0.04	-0.11 – 0.06	-0.53	0.595	0.03	0.07	-0.11 – 0.18	0.47	0.641
SN						0.05	0.05	-0.04 – 0.14	1.12	0.262
Session [Stress] * SN						-0.04	0.06	-0.16 – 0.08	-0.64	0.523
Run1 * SN						-0.02	0.04	-0.10 – 0.06	-0.44	0.660
Valence1 * SN						-0.02	0.04	-0.09 – 0.05	-0.46	0.644
(Session [Stress] * Run1) * SN						0.03	0.06	-0.08 – 0.15	0.57	0.569
(Session [Stress] Valence1) SN						-0.01	0.05	-0.12 – 0.09	-0.21	0.835
(Run1 * Valence1) * SN						0.02	0.04	-0.05 – 0.09	0.53	0.599
(Session [Stress] * Run1 * Valence1) * SN						-0.05	0.05	-0.15 – 0.06	-0.84	0.400
Random Effects
σ2	0.20					0.20
τ00	0.02 sub_nr					0.02 sub_nr
τ11	0.02 sub_nr.Session1					0.02 sub_nr.Session1
	0.00 sub_nr.Run1					0.00 sub_nr.Run1
ρ01	0.59					0.58
	-0.86					-0.92
N	63 sub_nr					63 sub_nr
Observations	444					444
Marginal R2 / Conditional R2	0.188 / NA					0.191 / NA

Plot: Main

plot.odd_valence <- ggplot(df_task.odd_val, aes(x=Scan, y=dprime, colour=Scan, fill=Scan)) + 
    geom_violin(alpha=0.5) + 
    geom_boxplot(width=0.1, alpha=0.2) +
    geom_hline(yintercept=0, colour="grey") +
    ggtheme2
plot.odd_valence + facet_grid(. ~ Valence)

Post-Hoc: Brain

joint_tests(models.odd_val[[2]], by=c("Session", "Run")) 

Session = Control, Run = Early:
 model term df1    df2 F.ratio p.value
 Valence      1 259.80   0.211  0.6464
 SN           1  99.13   0.295  0.5882
 Valence:SN   1 259.80   0.002  0.9638

Session = Stress, Run = Early:
 model term df1    df2 F.ratio p.value
 Valence      1 259.80   1.071  0.3017
 SN           1  94.03   0.198  0.6577
 Valence:SN   1 259.80   1.005  0.3171

Session = Control, Run = Late:
 model term df1    df2 F.ratio p.value
 Valence      1 259.80   1.298  0.2556
 SN           1 124.94   1.117  0.2926
 Valence:SN   1 259.80   0.499  0.4808

Session = Stress, Run = Late:
 model term df1    df2 F.ratio p.value
 Valence      1 259.80   0.075  0.7845
 SN           1  95.63   0.004  0.9497
 Valence:SN   1 259.80   0.001  0.9753

Hit+Miss Joint Plot

plot.odd_results <- ggarrange(
  ggarrange(ggarrange(NULL, (plot.odd_sess + xlab("") + labs(color="", fill="") + ylab("dPrime (a.u.)") + theme(axis.text.x=element_blank()) ),
                      NULL, (plot.odd_hits + xlab("") + ylab("Hits (count)")+theme(axis.text.x=element_blank()) ), 
                      NULL, (plot.odd_fas + xlab("") + ylab("False Alarms (count)")+theme(axis.text.x=element_blank()) ), 
                      ncol=6, common.legend = T, legend="bottom", widths=c(0.05, 1, 0.05, 1, 0.05, 1), labels =c("A", "", "B", "", "C")),
  (plot.sn_fas+ylab("False Alarms (count)")+xlab("SN Activity (a.u.)")), nrow=2, legend="bottom", labels=c("", "D")))

plot.odd_results
ggsave("figures/Figure_5_OddBeh.svg", device="svg", dpi =300, width=5, height = 6)

RT

df_task.oddrt <- df_task %>% subset(task==1 & valence==0 & rt!=0) %>% subset(rt>200) 
df_task.oddrt$valence <- (str_split_fixed(df_task.oddrt$stimulus, pattern="_", n =3))[,2]

Models

# Ma
signals=c("", "SN")
models.odd_rts <- mcmapply(signals, FUN=function(X){fx.par_model(signal=X,
                                               dv="rt",
                                               iv=ifelse(X=="",  "Session*Phase", "1+Session*Phase*"),
                                               re=ifelse(X=="",  "1+Session*Phase", paste("1+Session*Phase+", X, sep="")),
                                               sub_id="sub_id",
                                               fam=Gamma("log"),
                                               data=df_task.oddrt)}, mc.cores=5)

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be used

# Print HTML Output
asis_output(tab_model(models.odd_rts, dv.labels = c("Signal", "Signal ~ SN"), show.stat =T, show.se=T)$knitr)

	Signal					Signal ~ SN
Predictors	Estimates	std. Error	CI	Statistic	p	Estimates	std. Error	CI	Statistic	p
(Intercept)	603.46	5.97	591.87 – 615.27	647.37	<0.001	595.57	8.03	580.03 – 611.53	473.61	<0.001
Session [Stress]	0.97	0.01	0.96 – 0.99	-3.09	0.002	0.97	0.01	0.95 – 1.00	-1.65	0.098
Phase1	0.98	0.01	0.97 – 0.99	-4.09	<0.001	0.99	0.01	0.97 – 1.01	-1.22	0.223
Session [Stress] * Phase1	1.00	0.01	0.98 – 1.01	-0.69	0.492	0.99	0.01	0.96 – 1.01	-1.04	0.299
SN						1.01	0.01	0.99 – 1.03	1.29	0.196
Session [Stress] * SN						1.00	0.01	0.98 – 1.02	-0.27	0.787
Phase1 * SN						0.99	0.01	0.98 – 1.00	-1.41	0.160
(Session [Stress] Phase1) SN						1.01	0.01	0.99 – 1.03	0.75	0.455
Random Effects
σ2	0.06					0.06
τ00	0.01 sub_id					0.01 sub_id
τ11	0.00 sub_id.Session1					0.00 sub_id.Session1
	0.00 sub_id.Phase1					0.00 sub_id.Phase1
	0.00 sub_id.Session1:Phase1					0.00 sub_id.SN
						0.00 sub_id.Session1:Phase1
ρ01	-0.01					-0.01
	-0.09					-0.05
	-0.10					-0.58
						-0.26
ICC	0.08					0.08
N	70 sub_id					70 sub_id
Observations	11064					11064
Marginal R2 / Conditional R2	0.011 / 0.089					0.013 / 0.092

Post-Hoc: Main

emmeans(models.odd_rts[[1]],  list(pairwise ~Session*Phase), by=c("Phase"), adjust = "tukey")

$`emmeans of Session | Phase`
Phase = Early:
 Session emmean      SE  df asymp.LCL asymp.UCL
 Control  6.381 0.01084 Inf     6.360     6.402
 Stress   6.348 0.01101 Inf     6.327     6.370

Phase = Late:
 Session emmean      SE  df asymp.LCL asymp.UCL
 Control  6.424 0.01155 Inf     6.402     6.447
 Stress   6.401 0.01123 Inf     6.379     6.423

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$`pairwise differences of Session | Phase`
Phase = Early:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.0330 0.0114 Inf   2.895  0.0038

Phase = Late:
 2                estimate     SE  df z.ratio p.value
 Control - Stress   0.0233 0.0116 Inf   2.001  0.0454

Results are given on the log (not the response) scale. 

Plot: Main

df_task.oddrt_summ <- summarySEwithin(df_task.oddrt,idvar="sub_id", measurevar="rt", withinvars =c("Session", "Phase"))

Automatically converting the following non-factors to factors: Phase

df_task.oddrt_summ$Phase <- as.numeric(as.factor(df_task.oddrt_summ$Phase))
df_task.oddrt_summ$Session <- relevel(df_task.oddrt_summ$Session, "Stress")
ggplot(df_task.oddrt_summ, aes(x=Phase, y=rt, color=Session, fill=Session)) +
    geom_line() + 
    geom_errorbar(aes(ymin=rt-se, ymax=rt+se, color=Session),size=1, width=0.15, na.rm=T)+ 
  ggtheme2 + scale_x_continuous(breaks=c(1,2,3,4,5,6)) 

Post-Hoc: Brain

joint_tests(models.odd_rts[[2]], by=c("Session", "Phase"), adjust="none")

Session = Control, Phase = Early:
 model term df1 df2 F.ratio p.value
 SN           1 Inf   0.002  0.9634

Session = Stress, Phase = Early:
 model term df1 df2 F.ratio p.value
 SN           1 Inf   0.208  0.6484

Session = Control, Phase = Late:
 model term df1 df2 F.ratio p.value
 SN           1 Inf   3.336  0.0678

Session = Stress, Phase = Late:
 model term df1 df2 F.ratio p.value
 SN           1 Inf   0.967  0.3254

Mediation

library(mediation)

Loading required package: mvtnorm
Loading required package: sandwich
mediation: Causal Mediation Analysis
Version: 4.5.0


Attaching package: ‘mediation’

The following object is masked from ‘package:psych’:

    mediate

mediation.med <- lme4::lmer(SN ~ Session*Phase + (1|sub_nr), data=df_task.oddrt, REML=F)
mediation.main <- lme4::lmer(log10(rt) ~ Session*Phase + SN + (1|sub_nr), data=df_task.oddrt, REML=F)
mediation.analysis <- mediate(model.m=mediation.med, model.y=mediation.main, treat ="Session", control.value="Control", mediator = "SN", sims=5000)

treatment and control values do not match factor levels; using Control and Stress as control and treatment, respectively

summary(mediation.analysis)

Causal Mediation Analysis 

Quasi-Bayesian Confidence Intervals

Mediator Groups: sub_nr 

Outcome Groups: sub_nr 

Output Based on Overall Averages Across Groups 

                          Estimate 95% CI Lower 95% CI Upper p-value    
ACME (control)            3.53e-04     3.03e-05         0.00   0.029 *  
ACME (treated)            3.53e-04     3.03e-05         0.00   0.029 *  
ADE (control)            -1.01e-02    -1.38e-02        -0.01  <2e-16 ***
ADE (treated)            -1.01e-02    -1.38e-02        -0.01  <2e-16 ***
Total Effect             -9.75e-03    -1.35e-02        -0.01  <2e-16 ***
Prop. Mediated (control) -3.56e-02    -8.53e-02         0.00   0.029 *  
Prop. Mediated (treated) -3.56e-02    -8.53e-02         0.00   0.029 *  
ACME (average)            3.53e-04     3.03e-05         0.00   0.029 *  
ADE (average)            -1.01e-02    -1.38e-02        -0.01  <2e-16 ***
Prop. Mediated (average) -3.56e-02    -8.53e-02         0.00   0.029 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 11064 


Simulations: 5000 

FA by RT

Models

df_task.odd_rt_fa <- df_task.oddrt %>% dplyr::group_by(sub_nr, Session, Phase) %>% dplyr::summarise(rt=mean(rt)) %>% dplyr::rename(Run=Phase) %>% ungroup() %>% merge(., df_task.odd, by=c("sub_nr", "Session", "Run"))

`summarise()` has grouped output by 'sub_nr', 'Session'. You can override using the `.groups` argument.

model.odd_fa_rt <- lmer(rt ~ Session*Run*odd_FA + (1+Session+Run|sub_nr) , 
                        data=df_task.odd_rt_fa, 
                        contrast=list(Session="contr.treatment"), 
                        control=lmerControl(calc.derivs = F, optimizer="bobyqa", optCtrl=list(maxfun=100000)))
asis_output(tab_model(model.odd_fa_rt)$knitr)

	rt
Predictors	Estimates	CI	p
(Intercept)	631.60	596.33 – 666.86	<0.001
Session [Stress]	-33.31	-71.26 – 4.64	0.085
Run1	-43.64	-61.64 – -25.65	<0.001
odd_FA	-2.74	-6.10 – 0.62	0.110
Session [Stress] * Run1	29.11	2.47 – 55.75	0.032
Session [Stress] * odd_FA	2.20	-1.99 – 6.38	0.304
Run1 * odd_FA	3.16	0.93 – 5.40	0.005
(Session [Stress] * Run1) * odd_FA	-3.69	-6.69 – -0.69	0.016
Random Effects
σ2	817.56
τ00 sub_nr	8707.72
τ11 sub_nr.Session1	803.41
τ11 sub_nr.Run1	123.85
ρ01	0.11
	-0.66
ICC	0.92
N sub_nr	63
Observations	222
Marginal R2 / Conditional R2	0.049 / 0.920

check_model(model.odd_fa_rt)

`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

#Gaus 2670 #GausLog 3336.29 # Gamma 7770 #Gamma Log

Post-Hoc

joint_tests(model.odd_fa_rt, by=c("Session", "Run"), adjust="none")

Session = Control, Run = Early:
 model term df1    df2 F.ratio p.value
 odd_FA       1 109.60   0.033  0.8571

Session = Stress, Run = Early:
 model term df1    df2 F.ratio p.value
 odd_FA       1 107.68   0.355  0.5527

Session = Control, Run = Late:
 model term df1    df2 F.ratio p.value
 odd_FA       1 116.29   9.712  0.0023

Session = Stress, Run = Late:
 model term df1    df2 F.ratio p.value
 odd_FA       1 107.94   0.000  0.9934

Plot

emmip(model.odd_fa_rt, Session + Run ~ odd_FA, cov.reduce = range)+ ggtheme2# + scale_color_manual(values=stress_colors) 

rcorr(df_task.odd_rt_fa$rt, df_task.odd_rt_fa$odd_FA)

      x     y
x  1.00 -0.07
y -0.07  1.00

n= 222 


P
  x      y     
x        0.2951
y 0.2951       

plot.odd_fa_rt <- ggplot(df_task.odd_rt_fa, aes(x=odd_FA, y=rt)) + geom_smooth(method="lm" )
plot.odd_fa_rt + facet_grid(Session~Run)

5 System Info

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/cogaff/raytut/.conda/envs/r_env/lib/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mediation_4.5.0   sandwich_3.0-1    mvtnorm_1.1-2     doSNOW_1.0.19     snow_0.4-4        iterators_1.0.13  foreach_1.5.1     future_1.21.0     kableExtra_1.3.4 
[10] extrafont_0.17    stargazer_5.2.2   car_3.0-11        carData_3.0-4     sjPlot_2.8.10     ggpubr_0.4.0      emmeans_1.6.2-1   performance_0.8.0 jtools_2.1.3     
[19] optimx_2021-6.12  afex_1.0-1        lmerTest_3.1-3    lme4_1.1-27.1     Matrix_1.4-0      corrplot_0.90     ppcor_1.1         MASS_7.3-54       corrr_0.4.3      
[28] Hmisc_4.6-0       Formula_1.2-4     survival_3.2-13   lattice_0.20-45   psycho_0.6.1      psych_2.1.6       zoo_1.8-9         data.table_1.14.2 janitor_2.1.0    
[37] hablar_0.3.0      broom_0.7.9       plyr_1.8.6        forcats_0.5.1     stringr_1.4.0     dplyr_1.0.7       purrr_0.3.4       readr_2.1.1       tidyr_1.1.3      
[46] tibble_3.1.3      ggplot2_3.3.5     tidyverse_1.3.1   readxl_1.3.1      knitr_1.36       

loaded via a namespace (and not attached):
  [1] utf8_1.2.2          tidyselect_1.1.1    htmlwidgets_1.5.4   grid_3.6.1          lpSolve_5.6.15      munsell_0.5.0       codetools_0.2-18    effectsize_0.4.5   
  [9] withr_2.4.2         colorspace_2.0-2    highr_0.9           rstudioapi_0.13     robustbase_0.93-8   ggsignif_0.6.3      Rttf2pt1_1.3.9      listenv_0.8.0      
 [17] labeling_0.4.2      mnormt_2.0.2        bit64_4.0.5         farver_2.1.0        datawizard_0.2.1    parallelly_1.27.0   vctrs_0.3.8         generics_0.1.0     
 [25] xfun_0.25           R6_2.5.0            assertthat_0.2.1    scales_1.1.1        vroom_1.5.7         nnet_7.3-16         gtable_0.3.0        globals_0.14.0     
 [33] qqplotr_0.0.5       rlang_0.4.11        systemfonts_1.0.2   splines_3.6.1       rstatix_0.7.0       extrafontdb_1.0     checkmate_2.0.0     reshape2_1.4.4     
 [41] abind_1.4-5         modelr_0.1.8        backports_1.4.1     tools_3.6.1         ellipsis_0.3.2      RColorBrewer_1.1-2  ggridges_0.5.3      Rcpp_1.0.7         
 [49] base64enc_0.1-3     rpart_4.1-15        cowplot_1.1.1       haven_2.4.3         ggrepel_0.9.1       cluster_2.1.2       fs_1.5.0            magrittr_2.0.1     
 [57] openxlsx_4.2.4      reprex_2.0.1        tmvnsim_1.0-2       packrat_0.7.0       sjmisc_2.8.7        hms_1.1.0           evaluate_0.14       xtable_1.8-4       
 [65] pbkrtest_0.5.1      rio_0.5.27          sjstats_0.18.1      jpeg_0.1-9          gridExtra_2.3       ggeffects_1.1.1     compiler_3.6.1      crayon_1.4.1       
 [73] minqa_1.2.4         htmltools_0.5.2     mgcv_1.8-38         tzdb_0.1.2          lubridate_1.7.10    DBI_1.1.1           sjlabelled_1.1.8    dbplyr_2.1.1       
 [81] see_0.6.4           boot_1.3-28         cli_3.0.1           insight_0.14.4      pkgconfig_2.0.3     numDeriv_2016.8-1.1 foreign_0.8-71      xml2_1.3.2         
 [89] svglite_2.0.0       webshot_0.5.2       estimability_1.3    rvest_1.0.1         snakecase_0.11.0    digest_0.6.27       parameters_0.14.0   rmarkdown_2.11     
 [97] cellranger_1.1.0    htmlTable_2.3.0     curl_4.3.2          nloptr_1.2.2.2      lifecycle_1.0.0     nlme_3.1-153        jsonlite_1.7.2      viridisLite_0.4.0  
[105] fansi_0.5.0         pillar_1.6.2        fastmap_1.1.0       httr_1.4.2          DEoptimR_1.0-9      glue_1.4.2          bayestestR_0.10.5   zip_2.2.0          
[113] png_0.1-7           pander_0.6.4        bit_4.0.4           stringi_1.7.3       latticeExtra_0.6-29